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ABSTRACT 


A  program  was  undertaken  to  develop  a  system  of  models  for  describ¬ 
ing  transport  and  diffusion  of  smokes  and  agents  in  complex  terrain 
under  time-varying  meteorological  conditions.  The  ultimate  goal  is  to 
provide  the  probabilities  for  observing  concentrations  (point  or  path- 
integrated)  above  a  specified  threshold.  The  necessary  components  to 
achieve  this  goal  are: 


(1)  a  wind  model  for  describing  airflow  in  complex  terrain  using 
minimal  input  data 

(2)  a  transport  and  diffusion  model  capable  of  simulating  trans¬ 
port  and  diffusion  under  non-steady,  spatially-varying  condi¬ 
tions 


(3)  a  method  for  producing  realistic  small-scale  variabilities  in 
concentration  distributions 


(4)  a  method  for  estimating  the  required  probabilities  from  the 
small-scale  distributions. 


The  first  two  items  have  been  completed  and  are  described.  The  wind 
model  FORTRAN  source  code  is  included  in  an  appendix.  The  transport  and 
diffusion  model  was  previously  described,  and  a  Users’  Guide  has  been 
issued  separately. 


Lidar  cross-sections  through  smoke  plumes  were  analyzed  and  a 
preliminary  method  for  estimating  the  likelihood  of  finding  a  clear  path 
through  the  plume  is  discussed. 


The  literature  was  reviewed  to  find  a  framework  for  generating  the 
required  small-scale  concentration  distributions.  The  concept  of 
fractal  dimension  is  shown  to  have  considerable  promise.  An  extensive 
literature  review  and  bibliography  on  the  subject  of  scaling,  fractal 
dimension,  and  applications  to  atmospheric  processes  is  included  as  an 
appendix.  The  report  describes  the  research  necessary  to  apply  the 
identified  concepts  and  to  complete  the  desired  system  of  models. 
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I  INTRODUCTION 


This  final  report  under  ARO  Contract  No.  DAAG29-83-K-0009  summa¬ 
rizes  the  most  important  results  that  have  been  obtained  during  the 
course  of  this  project.  Detailed  information  is  provided  about  those 
research  results  that  have  not  been  described  previously  in  detail,  as 
was  done  in  the  publications  listed  in  Appendix  A.  This  is  in  accord¬ 
ance  with  the  Army  Research  Office  reporting  instructions,  which  also 
require  a  statement  of  the  problem  that  was  studied,  given  in  the  fol¬ 
lowing  section.  The  required  list  of  participants  in  the  project  is 
included  as  Appendix  B;  Appendices  C  and  D  present  research  accomplished 
on  this  project,  but  not  previously  reported  in  detail. 

Although  the  complex  terrain  airflow  model  has  been  described  in 
publ ications ,  no  guide  for  users  was  published.  Therefore,  Appendix  C 
has  been  included  to  provide  an  abbreviated  set  of  instructions,  and  a 
listing  of  the  source  code  for  the  wind  model.  Similarly,  the  review  of 
the  applications  of  fractals  to  the  study  of  atmospheric  processes,  and 
the  bibliography  on  the  same  topic  that  was  compiled  have  not  been 
published  to  date.  For  completeness,  they  appear  here  as  Appendix  D. 


II  STATEMENT  OF  THE  PROBLEM  STUDIED 


This  study  has  addressed  the  problem  of  modeling  the  behavior  of 
smoke  plumes  or  agents  in  complex  terrain  during  temporally  and  spa¬ 
tially  varying  atmospheric  conditions.  In  addition  to  the  most  probable 
concentration  distributions ,  an  analytical  framework  was  sought  which 
could  be  used  to  estimate  the  probabilities  of  encountering  concentra¬ 
tions  (point  or  path-integrated)  exceeding  some  specified  value.  This 
is  a  large  and  difficult  problem,  and  was  divided  into  two  major  parts. 

The  first  part  of  the  problem  was  to  develop  a  deterministic  trans¬ 
port  and  diffusion  model  that  could  be  used  to  describe  the  "average"  or 
"most  probable"  concentration  distribution.  A  non-steady-state  Lagran- 
gian  model  was  developed  that  is  very  efficient  and  practical  (Ludwig, 
198^3,6,  1985a).  However,  it  became  apparent  during  the  course  of  the 
project  that  transport  and  diffusion  models  were  no  better  than  the  wind 
field  and  meteorological  models  used  to  provide  inputs.  Therefore,  a 
wind  model  was  developed  to  use  minimal  observational  data,  and  to  be 
efficient  and  practical  enough  for  field  use  (Ludwig  et  al ,  1985). 

The  more  difficult  part  of  the  problem  was  to  identify  a  framework 
that  could  be  used  to  estimate  concentration  probabilities.  Two  comple¬ 
mentary  task  efforts  were  undertaken  to  attack  this  problem.  First, 
relevant  data  from  past  field  tests  and  experiments  were  reviewed  and 
analyzed  in  an  attempt  to  develop  a  statistical  framework  that  could  be 
applied  to  real  data  to  define  the  probability  distributions  that  were 
sought.  The  existing  data  were  not  completely  adequate  for  this  pur¬ 
pose,  but  they  were  useful  for  developing  promising  approaches  and 
determining  what  further  information  might  be  needed.  This  approach 
presumed  that  the  spatial  distribution  of  material  could  be  known  in 
considerable  detail,  a  presumption  that  cannot  in  fact  be  met,  except  in 
the  statistical  sense.  This  leads  to  the  second  task  where  we  sought 
methods  to  describe  small  scale  variabilities  and  concentrations  realis¬ 
tically  and  in  a  fashion  consistent  with  the  larger  scale  features 
determined  by  the  genera]  airflow  and  meteorological  conditions. 

In  an  effort  to  develop  methods  for  defining  relatively  small-scale 
variabilities,  the  literature  was  reviewed,  and  a  promising  approach  was 
identified.  This  approach  is  based  on  the  use  of  fractals  (e.  g. 
Lovejoy  and  Schertzer,  1986).  There  is  evidence  that  fractals  can  be 
used  to  provide  the  bridge  between  large-scale  atmospheric  features  that 
can  be  modeled  deterministically  and  the  important  statistical  proper¬ 
ties  of  the  small-scale  features  that  cannot.  Although  we  have  not  been 
able  to  develop  the  complete  statistical  model,  the  information  neces¬ 
sary  to  complete  the  work  is  clearly  defined. 

In  summary,  the  complete  system  that  we  are  seeking  will  include 
the  following  components: 


(1)  A  wind  field  generator  that  can  be  used  with  limited  inputs 
for  generating  3-dimensional  winds  in  complex  terrain 

(2)  A  deterministic  transport  and  diffusion  model  that  can  define 
the  general,  most  probable  concentration  distribution  associ¬ 
ated  with  the  plume  or  cloud  of  material 

(3)  A  statistical  model  for  generating  realistic  small-scale 
concentration  fluctuations  within  the  plume 

(4)  A  statistical  model  for  calculating  probability  distributions 
for  point  or  path-integrated  concentrations  given  the  small- 
scale  concentration  pattern. 

Models  were  developed  on  this  project  that  meet  the  requirements  of  the 
first  two  items  in  the  above  list.  The  problems  associated  with  (4) 
were  not  completely  solved,  but  enough  progress  was  made  so  that  it  is 
believed  that  there  will  be  no  fundamental  difficulties  in  arriving  at 
such  a  model.  The  third  item  has  proven  to  be  by  far  the  most  diffi¬ 
cult,  but  we  feel  that  an  approach  involving  the  application  of  fractals 
can  be  successful  in  providing  the  necessary  information. 

The  remainder  of  this  report  discusses  significant  research 
results.  The  findings  with  regard  to  the  first  two  items  above  are  more 
complete,  and  papers  and  reports  have  already  been  prepared  describing 
results.  The  other,  more  preliminary  findings,  are  described  in  some¬ 
what  more  detail  because  papers  have  not  yet  been  prepared.  Although 
these  findings  are  still  preliminary  and  tentative,  we  believe  them  to 
be  quite  significant. 


Ill  SUMMARY  OF  IMPORTANT  RESULTS 


A.  Wind  Model  Description 


As  noted  before,  transport  and  diffusion  models  require  accurate 
specification  of  the  wind  fields.  However,  practical  applications, 
especially  those  of  greatest  interest  to  the  Army,  will  generally  be 
restricted  to  the  use  of  data  from  widely  separated  sites.  Furthermore, 
applications  in  the  field  must  usually  make  do  with  very  restrictive 
computational  facilities.  The  model  described  in  Appendix  C  was  devel¬ 
oped  to  provide  realistic  wind  fields  using  minimal  input  data,  and 
calculations  that  do  not  require  a  mainframe  computer. 

The  model  is  described  completely  (including  a  F0RTRAN77  source 
code)  in  Appendix  C.  Very  briefly,  we  attempted  to  devise  a  scheme  that 
made  maximum  use  of  the  available  observations  and  theoretical  and 
empirical  relationships  governing  atmospheric  processes.  The  result  is 
essentially  a  wind  interpolation  scheme  that  is  constrained  to  satisfy 
the  principal  of  conservation  of  mass  (nondi vergence) ,  and  an  empirical 
relationship  based  on  conservation  of  energy.  The  latter  constraint 
arises  from  the  fact  that  there  is  a  buoyant  restoring  force  in  a  stably 
stratified  atmosphere  whenever  air  is  displaced  from  its  equilibrium 
position.  Air  displaced  upward  will  be  cooler  (denser)  than  its  sur¬ 
roundings  and  subjected  to  a  downward  force,  while  downward-moving  air 
will  have  an  upward  restoring  force.  Work  is  required  to  move  a  volume 
of  air  from  its  equilibrium  altitude,  and  this  results  in  an  increase  in 
the  potential  energy  of  the  displaced  air.  Hunt  and  Snyder  (1980) 
reasoned  that  the  potential  energy  was  gained  at  the  expense  of  the 
kinetic  energy  so  that  for  a  given  wind  speed  and  temperature  stratifi¬ 
cation,  there  is  some  maximum  upward  displacement  that  can  be  induced  by 
the  underlying  terrain.  There  will  be  a  streamline  for  which  this 
upward  displacement  causes  the  air  to  rise  just  to  the  height  of  a 
terrain  obstacle.  This  is  the  "critical  streamline."  Air  below  this 
altitude  passes  around  the  obstacle.  Air  above  the  critical  streamline 
flows  over  the  obstacle.  The  validity  of  the  critical  streamline  con¬ 
cept  seems  to  be  well  supported  by  field  experiments  (Lavery  et  al, 
1982;  Egan,  1984). 

In  its  original  form,  the  critical  streamline  was  used  to  determine 
whether  or  not  a  plume  centerline  would  pass  over  a  hill  or  around  it. 
The  model  described  in  Appendix  C  has  extended  the  concept  to  the  defi¬ 
nition  of  three-dimensional  flow  fields.  It  is  used  to  define  the  shape 
of  "flow  surfaces"  within  which  the  streamlines  lie.  These  surfaces  may 
intersect  a  terrain  feature  or  pass  over  it.  If  they  intersect,  the 
principle  of  mass  conservation  (nondivergence)  forces  the  airflow  to 
pass  around  the  obstacle.  This  relatively  simple  concept  has  been 
implemented  in  the  model  described  in  Appendix  C. 

The  model  achieves  the  objectives  of  using  minimal  data.  It  uses  a 
vertical  profile  of  temperature.  It  is  likely  that  data  from  my  single 
Doppler  acoustic  sounder  could  be  substituted  for  the  temperature 


profile,  presuming  that  the  vertical  fluctuation  in  the  wind  can  be 
related  to  the  lapse  rate.  The  model  could  be  run  with  only  a  single 
wind  input.  However,  much  better  results  will  be  obtained  if  three  or 
more  surface  observations  are  available  so  that  estimates  of  upper-level 
winds  can  be  obtained  from  the  geostrophic  and  thermal  wind  relation¬ 
ships.  Another  alternative  would  be  to  have  Doppler  acoustic  sounder, 
pilot-balloon,  or  other  estimates  of  winds  aloft.  In  summary,  the  model 
requires  some  estimate  of  stability  stratification  of  the  atmosphere  and 
can  use  whatever  wind  observations  are  available.  Obviously,  its  accui — 
acy  increases  with  more  input  data. 

The  objective  of  computational  efficiency  has  been  met.  The 
samples  given  in  Appendix  C  required  only  about  15  seconds  of  central 
processing  unit  (CPU)  time  on  a  VAX  11/782  computer.  The  practicality 
of  the  model  is  demonstrated  by  the  fact  that  it  has  been  incorporated 
into  a  larger  code  being  developed  at  the  U.S.  Army’s  Atmospheric 
Sciences  Laboratory  at  White  Sands,  New  Mexico  (Dr.  R.  Meyers,  personal 
communication,  1986).  It  is  understood  that  this  latter  application  is 
on  an  IBM- AT  computer. 


B.  Transport  and  Diffusion  Model 

The  transport  and  diffusion  models  developed  on  this  project  have 
been  described  elsewhere.  Two  versions  have  been  written.  The  first 
version  was  written  in  BASIC  for  use  on  an  Apple-II  or  IBM-PC  (Ludwig, 
1 98^ ) .  That  version  was  upgraded  and  converted  to  a  second  version  in 
FORTRAN77 .  A  complete  users’  guide  was  prepared  for  the  FORTRAN77 
version  as  a  Technical  Project  Report  (Ludwig,  1985). 

The  transport  and  diffusion  models  were  developed  before  the  wind 
model  described  in  the  preceding  section.  As  a  result,  they  were 
designed  to  use  the  outputs  from  models  that  invoked  only  the  mass- 
consistency  constraint.  Small  computers  were  slower  and  had  more  severe 
memory  constraints  at  the  time  the  models  were  developed  than  they  do 
now,  so  attempts  were  made  to  reduce  the  computations  required  for  the 
wind  fields.  The  method  described  by  Ludwig  and  Byrd  (1980)  for 
exploiting  the  linear  properties  of  mass-consistent  wind  models  was 
employed  so  that  most  of  the  calculations  could  be  done  ahead  of  time  on 
a  large  machine,  and  then  used  to  generate  the  winds  from  linear  com¬ 
binations  of  those  precalculated  solutions.  The  speed  and  available 
memory  of  small  computers  is  now  such  that  the  transport  and  diffusion 
models  developed  earlier  should  probably  be  modified  so  that  the  winds 
are  calculated  online,  using  the  model  described  in  the  preceding  sec¬ 
tion.  Time  was  not  available  to  make  this  modification,  although  it 
should  not  be  particularly  difficult. 

The  models  use  a  Lagrangian  puff  approach  to  describe  a  continuous 
emission.  Each  puff  has  a  bivariate  Gaussian  concentration  distribu¬ 
tion.  The  puffs  are  moved  according  to  the  winds  at  their  centers. 
They  expand  according  to  the  atmospheric  stability  at  the  time  they  are 


advected.  As  expansion  causes  the  puffs  to  overlap,  they  are  combined 
to  reduce  required  computation. 

Concentrations  over  a  grid  of  points  at  the  surface  are  calculated 
hourly  by  summing  the  contributions  from  all  puffs  within  about  three 
standard  deviations  of  the  point.  The  usual  assumptions  of  "reflection" 
at  the  surface  are  used. 

Details  of  the  theory  and  its  implementation  can  be  found  in  two 
publications  prepared  during  the  course  of  this  contract  (Ludwig,  198-4, 
1985). 


C .  Probability  Analysis  Framework  for  Point  and  Path  Integrated 

Concentrations 

The  framework  that  has  evolved  for  assessing  probabilities  of 
encountering  various  concentration  levels  (or  path- integrated  concentra¬ 
tions)  is  a  three  step  process: 

(1)  Average  (or  most  probable)  concentration  patterns  are  calcu¬ 
lated  from  the  deterministic  wind  and  dispersion  models 
described  above 

(2)  A  subgrid  scale  spatial  concentration  pattern  is  superimposed 
on  the  smooth  pattern  from  Step  1 

(3)  Point  or  line-of-sight  concentration  probabilities  are  esti¬ 
mated  either  from  repeated  realizations  of  Step  2,  or  from 
values  at  similar  locations  or  along  similar  paths  from  a 
single  subgrid-scale  pattern. 

The  generation  of  realistic  concentration  patterns  appears  possible 
by  adapting  existing  methods  for  generating  fields  of  specified  fractal 
dimension  (Pentland,  1984;  Lovejoy  and  Mandelbrot,  1985;  Medler,  Gelberg 
and  Burkhart,  1986;  Wilson  et  al,  1986).  As  noted  by  Medler  et  al 
(1986),  this  requires  knowledge  of  the  fractal  dimension  of  the  patterns 
of  smoke  concentration  in  a  plume.  Securing  this  missing  information 
will  be  a  major  research  effort.  Section  III-D  describes  the  required 
research  program  that  was  developed  during  this  contract. 

The  relevance  of  fractals  to  this  problem  is  described  in  Appendix 
D  which  reviews  recent  developments  in  the  use  of  fractals  to  describe 
spatially  inhomogeneous  scalar  and  vector  atmospheric  fields.  Several 
methods  are  described  for  estimating  fractal  dimensions  and  generating 
the  corresponding  fields.  Since  that  review  was  written,  other  methods 
have  been  identified  which  seem  even  more  suitable  for  the  anisotropic 
atmospheric  distributions  (e.g.  Wilson  et  al,  1986;  Schertzer  and 
Lovejoy,  1986). 

Figure  1  shows  an  example  of  one  approach  that  has  been  taken  to 
estimate  clear  lines-of-sight  from  observed  lidar  cross-section  data. 
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Points  have  been  plotted  showing  the  distance  from  the  center  of  mass  of 
a  smoke  plume  on  the  abcissa  versus  distance  to  the  nearest  pixel  whose 
value  (backscatter )  is  above  a  selected  value.  In  Figure  1(a),  only 
lines-of-sight  that  are  within  30°  of  the  direction  toward  the  center  of 
the  plume  have  been  included.  The  picture  elements  are  nearly  square, 
about  20  m  on  a  side.  The  clustering  of  points  at  the  top  of  the  plot 
are  those  that  pass  unobstructed  beyond  the  domain.  As  expected,  the 
open  paths  are  shorter  near  the  center  of  the  plume.  Figures  1(b), 
1(c),  and  1(d)  show  that  it  is  increasingly  likely  that  a  clear  line-of- 
sight  will  be  found  in  directions  away  from  the  plume  center.  The  data 
in  Figure  1  could  be  converted  easily  to  a  probability  function.  Back¬ 
scatter  values  along  the  paths  could  also  be  summed  to  obtain  path- 
integrated  values  whose  probabilities  of  occurence  could  be  estimated. 


D.  Description  of  Experimental  and  Theoretical  Research  Required  to 
Implement  a  Fractal-Based  Probabilistic  Methodology 


A  research  plan  was  developed  during  the  course  of  this  contract  to 
provide  a  description  of  the  effort  needed  to  implement  the  probabilis¬ 
tic  framework  described  above.  Four  main  tasks  are  required: 

(1)  Acquisition  of  data 

(2)  Development  or  refinement  of  fractal  analysis  techniques 

(3)  Data  analysis  and  interpretation  using  fractal  methodology 

(4)  Adaptation  of  available  techniques  or  development  of  new  ones 
to  generate  fields  with  defined  fractal  and  larger-scale 
characteristics. 

These  are  discussed  in  more  detail  in  the  following  sections. 


1 .  Data  Acquisition 

Analysis  techniques  for  digitized  lidar  cross-sections  have  been 
developed  during  the  present  effort.  These  techniques  make  use  of  a  PDP 
11/23  computer  with  a  Peritech  graphics  package  that  facilitates  visual 
interpretation  of  the  results.  The  digitized  lidar  backscatter  cross- 
sections  have  proven  to  be  very  amenable  to  automated  analysis,  but  the 
currently  available  data  have  some  shortcomings,  namely: 

•  They  have  limited  spatial  resolution,  about  15  to  30  m  horizon¬ 
tal  by  3  to  6  m  vertical 

•  Most  of  the  available  data  are  for  elevated,  buoyant  plumes  and 
far  downwind  distances. 

The  first  shortcoming  is  not  fundamental  and  could  be  corrected  by 
operating  the  equipment  differently.  It  should  be  possible  to  achieve  a 
resolution  of  about  3  x  3  m. 
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The  second  shortcoming  requires  appropriate  field  tests  which  meet 
the  following  criteria: 

•  Ground  level  release  of  smoke  or  other  aerosol  tracer  material 
(not  dense  enough  to  be  opaque  to  the  lidar) 

«  Higher  resolution  backscatter  cross  sections  measured  at  2  or  3 
distances  closer  (about  0.5  to  2.5  km)  to  the  source  than  cur¬ 
rently  available  data 

•  Several  meteorological  conditions  represented,  e.g.  stable 
mornings,  convective  afternoons,  neutral,  transition  periods, 
etc. 

Airborne  lidar  systems  are  available  that  could  collect  digitized 
data  with  a  resolution  of  a  few  meters  as  required  (Uthe  et  al,  1980). 

The  availability  and  suitability  of  data  from  the  stereo,  multi- 
spectral  imaging  systems  (formerly  referred  to  as  MIDAS)  that  are  oper¬ 
ated  at  some  Smoke  Week  tests  (e.g.  Farmer  et  al,  1986)  should  also  be 
explored.  If  these  data  could  be  used,  they  might  obviate  the  need  for 
the  special  tests  described  above,  especially  if  it  is  possible  to 
select  planes  through  the  smoke  of  different  orientation  so  that  the 
degree  of  isotropy  of  the  patterns  could  be  tested  and  some  of  the  newer 
fractal  analysis  techniques  (Schertzer  and  Lovejoy,  1986)  could  be 
applied. 

The  output  velocity  fields  from  large  eddy  computer  simulations  are 
another  potential  source  of  valuable  information.  However,  there  are 
some  significant  difficulties  involved  with  using  these  data.  Foremost 
among  these  is  the  problem  of  storing  and  retrieving  the  massive  volume 
of  results  that  are  generated,  i.e.  literally  millions  of  numbers  for 
every  few  seconds  of  simulated  time.  Although  this  approach  3eems 
attractive  because  of  the  ability  to  generate  true  ensembles  of  cases, 
it  is  not  clear  that  the  problems  are  surmountable.  If  the  data  were 
usable,  "experiments"  could  be  conducted  to  cover  a  variety  of  condi¬ 
tions,  e.g.: 

•  Different  sizes  and  shapes  of  initial  particle  clouds 

•  Continuous  releases  in  a  field  with  mean  transport 

•  Different  locations  if  an  inhomogeneous  field  is  used. 

Such  computer  experiments  have  the  potential  to  identify  relation¬ 
ships  between  observable  patterns  of  smoke  diffusion  and  the  charac¬ 
teristics  of  atmospheric  turbulence.  Chorin's  (1982)  analysis  of  the 
evolution  of  the  vorticity  field  associated  with  a  bent  vortex  tube 
provides  considerable  encouragement  for  this  kind  of  analysis. 


2.  Develop  Fractal  Analysis  Methods  for  _ Inhomogeneous _ and 

Anisotropic  Fields  ’ 

Only  two  of  several  available  fractal  estimation  methods  have  been 
developed  and  used  during  the  current  research  effort,  and  those  two 
methods  were  used  under  the  assumption  that  the  processes  were  iso¬ 
tropic.  Lovejoy  and  Schertzer  (1986)  have  recently  pointed  out  the 
importance  of  anisotropy  in  the  atmosphere.  Some  of  the  methods  avail¬ 
able  for  estimating  fractal  dimensions  are: 

(1)  A  subcell  counting  method  (e.g.  Mandelbrot,  1975;  Ludwig  and 
Nitz,  1986) 

(2)  An  average  power  spectrum  method  (e.g.  Pentland,  1 984 ) 

(3)  A  method  based  on  the  probability  distribution  of  the  differ¬ 
ences  in  the  variable  versus  separation  distance  (e.g. 
Schertzer  and  Lovejoy,  1983) 

(M)  A  method  based  on  the  relationship  between  the  area  and  the 
perimeter  of  regions  enclosed  by  isopleths  (Lovejoy  and 
Mandelbrot,  1985) 

(5)  A  method  for  estimating  "elliptical  fractal  dimension" 
(Lovejoy  and  Schertzer,  1986)  using  a  variation  of  the  subcell 
counting  method  (Schertzer  and  Lovejoy,  1986). 

This  is  a  rapidly  evolving  area  of  research  and  it  will  be  essen¬ 
tial  that  new  methods  be  examined  for  applicability  as  they  are  devel¬ 
oped.  The  fifth  method  listed  above  shows  considerable  promise  because 
of  its  ability  to  describe  anisotropic  features.  It  appears  to  be 
particularly  applicable  to  3-d  data  like  those  from  the  stereo  imaging 
systems . 


3.  Data  Analysis  and  Interpretation  with  Fractal  Methods 

Although  several  of  the  methods  discussed  above  should  be  mathema¬ 
tically  equivalent,  there  are  some  differences  that  arise  from  the  way 
that  they  are  applied,  so  a  comparison  of  results  is  warranted.  The 
most  important  of  the  differences  is  between  the  cell  counting  method 
and  the  others.  In  the  isotropic  case,  the  cell  counting  method  defines 
a  fractal  dimension  for  a  topologically  one-dimensional  surface,  the 
isopleth.  It  is  thought  that  the  intersection  of  a  surface  and  a  plane 
has  a  fractal  dimension  one  less  than  that  of  the  surface  itself 
(Mandelbrot,  1975)  so  it  would  be  easy  to  reconcile  the  two  approaches 
except  that  isopleths  tend  to  be  arranged  in  concentric  patterns  about 
the  plume  center.  Therefore,  the  fractal  dimensions  of  high-valued 
isopleths  will  be  more  repr esentati ve  of  the  interior  of  the  plume  than 
are  the  fractal  dimensions  of  the  lower  value  isopleths.  The  other 
methods  will  have  to  be  applied  to  fairly  large  subregions  of  the  plume 
in  order  to  have  samples  that  are  sufficiently  large  to  characterize  the 


probability  functions  or  spectral  densities.  These  subregions  are 
likely  to  cut  across  a  range  of  isopleth  "rings". 

The  discussion  of  possible  methodological  differences  in  the  pre¬ 
ceding  section  points  out  the  potential  for  differences  from  one  part  of 
a  plume  to  another.  Without  systematic  examination,  it  is  not  clear 
whether  or  not  these  differences  will  be  found.  At  first,  it  seems  as 
though  clean  air  entrainment  at  the  edge  of  a  plume  might  lead  to  a 
different  pattern  than  in  the  interior.  However,  this  may  only  cause 
differences  in  the  amplitude  of  fluctuations,  but  their  spatial  pattern 
may  be  the  same.  Stated  differently,  the  slope  of  the  spectral  density 
function  may  be  the  same  although  the  amplitudes  may  be  larger  at  the 
edge  than  in  the  middle. 

During  periods  of  fumigation  or  lofting,  part  of  the  plume  will  be 
in  a  stable  layer  and  part  in  a  more  turbulent  zone.  Differences, 
especially  in  the  nature  of  any  anisotropic  effects  are  quite  likely  and 
should  be  investigated. 

The  overall  dimension  of  the  plume  almost  certainly  will  set  some 
upper  limit  on  the  range  over  which  scaling  occurs.  As  noted  earlier, 
there  are  apt  to  be  differences  between  the  vertical  and  the  horizontal. 
The  relationships  between  scaling  range  and  plume  dimension  should  be 
examined. 

Lidar  cross  sections  taken  through  the  widely  dispersed  plumes  in 
convective  boundary  layers  have  a  very  different  appearance  from  the 
more  compact  plumes  found  when  there  is  less  convective  activity.  It  is 
not  certain  whether  or  not  these  cases  have  a  different  fractal  dimen¬ 
sion,  but  at  the  very  least,  it  seems  reasonable  to  assume  that  the 
"spheroscale"*  will  be  changed  by  the  presence  or  absence  of  convection. 

If  one  assumes  that  the  eddy  scales  affect  the  patterns  of  smoke 
distribution  within  a  plume,  then  there  is  reason  to  believe  that  there 
will  be  variations  with  vertical  position  in  the  mixed  layer.  Near  the 
upper  and  lower  boundaries  of  the  convective  mixing  layer,  vertical 
motions  are  inhibited,  so  local  patterns  may  be  vertically  compressed. 
One  obvious  objective  of  an  analysis  of  the  vertical  changes  in  fractal 
dimension  within  the  mixing  layer  would  be  to  see  if  the  spheroscale  is 
a  function  of  the  mixing  height.  An  even  more  important  objective  would 
be  to  relate  the  properties  of  smoke  patterns  to  the  statistics  of  the 
eddy  motions  from  large  eddy  simulation  and  field  experiments. 


"Spheroscale"  is  the  term  used  by  Lovejoy  and  Schertzer  (1986)  for  the 
scale  at  which  the  characteristic  vertical  and  horizontal  dimensions  of 
the  pattern  are  equal. 


From  the  practical  standpoint,  an  ability  to  relate  the  fractal 
properties  of  smoke  distributions  to  the  prevailing  meteorology  would  be 
most  important.  It  is  doubtful  that  this  could  be  done  for  the  complete 
range  of  meteorological  conditions  because  of  limitations  on  the  condi¬ 
tions  under  which  lidar  data  can  be  collected.  However,  there  already 
exists  a  body  of  data  covering  a  broad  range  of  atmospheric  stability, 
wind  speed  and  surface  heating  conditions.  Future  experiments  should 
extend  these  to  other  conditions  and  various  positions  within  the  mixed 
layer. 

If  recent  studies  of  turbulent  diffusion  modeling  can  serve  as  a 
guide,  then  mixing  depth  (H),  Monin-Obhukov  length  (L),  convective 
scaling  velocity  (w#),  friction  velocity  (u*),  wind  fluctuations 
(Og  and  o  )  and  eddy  heat  flux  are  good  candidate  parameters  for  charac¬ 
terizing  \he  meteorological  conditions.  For  this  reason,  it  would  be 
preferable  to  use  data  sets  for  which  the  appropriate  special  corollary 
meteorological  measurements  were  available.  However,  such  is  not  always 
the  case,  so  it  will  often  be  necessary  to  estimate  the  appropriate 
parameters  from  more  conventional  observations  by  following  quidelines 
such  as  those  suggested  by  Holtslag  et  al  (1985)  or  Hanna  (19814). 


4 .  Develop  Appropriate  Fractal  Simulation  Methods 

The  next  step  in  the  process  will  be  to  take  the  relationships 
gleaned  from  the  analysis  described  in  the  preceding  section  and  use  it 
as  the  basis  for  specific  fractal  generation  algorithms.  It  is  not 
necessary  to  develop  entirely  new  approaches,  but  rather  to  use  the 
fractal  characteristics  derived  from  analysis  and  interpretation  with 
existing  methods  (e.g.  Lovejoy  and  Mandelbrot,  1985;  Norton,  1982;  Voss, 
1982;  Wilson  et  al ,  1986).  It  will  be  necessary  to  superimpose  the 
fractal  patterns  on  the  general  distributions  of  concentration,  such  as 
the  frequently  used  Gaussian  distribution. 

A  problem  may  be  encountered  if  fractal  characteristics  vary  over 
the  area  of  a  plume,  as  preliminary  analysis  suggests  may  be  the  case. 
Most  existing  image  generation  algorithms  have  only  to  deal  with  pat¬ 
terns  and  textures  that  have  reasonably  well  defined  edges  rather  than 
smooth  transitions  from  one  to  another.  The  random  cylinder  method  of 
Lovejoy  and  Mandelbrot  (1985)  appears  promising  if  modified  to  use 
elliptical  cross-sections  to  account  for  anisotropy.  The  parameters  of 
the  random  generating  function  would  be  functions  of  location  allowing 
for  the  desired  smooth  transitions. 

Once  the  appropriate  spatial  patterns  can  be  generated,  it  will  be 
possible  to  combine  them  with  the  results  from  more  conventional  models 
(e.g.  complex  terrain  airflow  and  nonsteady  state  transport  and  diffu¬ 
sion  models)  to  generate  more  realistic  spatial  distributions  of  concen¬ 
tration.  These  distributions  can,  in  turn,  serve  as  a  basis  for  esti¬ 
mates  of  the  probability  of  exceeding  specified  point  or  path  integrated 
concentrations .  Under  the  assumption  that  patterns  remain  "frozen"  for 
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I  INTRODUCTION 


As  a  general  rule,  sites  from  which  wind  observations  are  available 
are  widely  separated.  If  the  atmosphere  is  reasonably  well  mixed  and 
the  terrain  is  flat,  simple  interpolation  techniques  can  be  used  to 
estimate  winds  at  any  given  location  between  the  observing  sites.  There 
are  many  applications  for  which  wind  fields  obtained  by  simple  interpo¬ 
lation  techniques  are  not  adequate.  Any  application  that  involves  the 
calculation  of  fluxes  and  concentrations  will  suffer  if  the  wind  field 
has  appreciable  divergence.  The  most  common  of  such  applications  is  in 
the  modeling  of  transport  and  diffusion.  A  divergent  wind  field  can 
completely  obliterate  the  effects  of  turbulent  diffusion,  chemical 
reactions,  and  even  transport  in  such  models. 

When  complex  terrain  and  stable  atmospheric  stratifications  are 
introduced,  simple  interpolation  schemes  can  be  quite  misleading.  Some 
improvements  can  be  achieved  by  invoking  the  principle  of  conservation 
of  mass  to  produce  nondivergent  wind  fields,  but  there  still  remains  the 
question  of  whether  air  will  flow  over,  or  around  an  obstacle.  The 
result  obtained  by  a  transport  and  diffusion  model  is  very  dependent 
upon  the  choice  that  is  made.  The  effects  extend  beyond  determining 
which  areas  are  affected  by  a  plume;  the  magnitude  of  the  effects  is 
also  determined  by  the  wind  field.  For  example,  a  plume  that  impacts  a 
hill  will  produce  much  greater  ground  level  concentrations  than  one  that 
passes  to  one  side  and  stays  well  above  the  surface.  There  is  an 
obvious  need  for  wind  interpolation  schemes  that  can  estimate  winds  on  a 
three-dimensional  grid,  especially  objective  methods  with  few  "adjust¬ 
able  constants"  and  modest  computational  requirements . 

An  ideal  approach  would  invoke  all  the  laws  of  thermodynamics  and 
fluid  motion  in  the  interpolation  process.  In  essence,  this  would  be  a 
planetary  boundary  layer  model  that  incorporated  all  the  available 
observations.  It  is  likely  that  it  would  be  at  least  as  complicated  as 
existing  planetary  boundary  layer  models,  so  that  it  could  not  be  widely 
used  with  inexpensive  machines.  This  paper  limits  itself  to  a  more 
modest  approach.  The  principle  of  conservation  of  mass  (nondivergence) 
has  been  invoked  in  several  wind  interpolation  schemes  (e.g.,  Sherman, 
1978;  Dickerson,  1978;  Bhumralkar  et  al.,  1980;  Goodin  et  al.,  1980;  and 
Endlich  et  al  ,  1982).  The  next  logical  step  seems  to  be  the  incorpora¬ 
tion  of  conservation  of  energy  to  account  for  the  changes  in  potential 
energy  with  vertical  motion  when  the  atmosphere  is  stably  stratified. 

Current  approaches  to  wind  field  interpolation  involve  two  steps. 
First,  winds  are  estimated  at  the  points  in  a  2-  or  3“dimensional  grid 
by  some  standard,  unconstrained  interpolation  technique.  The  second 
step  is  to  adjust  this  initial  wind  field  so  that  it  satisfies  some 
specified  constraints.  The  initial  wind  field  is  typically  obtained  by 
linear  interpolation,  inverse  distance  (or  inverse  distance  squared) 
weighting  or  a  Gaussian  weighting  scheme.  Most  commonly,  the  con¬ 
straints  involve  conservation  of  mass.  Sherman  (1978)  and  Bhumralkar  et 
al  (1980)  both  further  constrain  the  adjustment  so  that  it  represents  a 


"minimum"  change  in  the  initia]  field  required  to  reduce  the  divergence 
(and  convergence)  below  some  specified  level.  Endlich  (1967)  has  a  more 
general  technique  where  the  initial  wind  fields  can  be  adjusted  to 
satisfy  seme  specified  fields  of  divergence  and  vertical  vorticity 
component.  These  approaches  provide  a  useful  starting  point.  Effects 
of  stable  stratification  on  air  flow  are  discussed  in  the  next  section 
to  provide  the  rationale  for  the  method  described  later. 


II  BUOYANCY  EFFECTS 


The  following  discussion  deals  with  those  cases  where  potential 
temperature  increases  with  height,  i.e.,  30/3z>O.  Furthermore,  we  shall 
only  deal  with  the  situations  where  the  atmospheric  processes  are 
approximately  adiabatic.  Thus,  the  following  discussion  assumes  that  no 
condensation  or  evaporation  takes  place  and  the  time  scales  are  of  the 
order  of  a  few  hours  or  less,  so  that  there  will  be  a  buoyant  restoring 
force  whenever  the  air  is  displaced  from  its  equilibrium  position.  Air 
displaced  upward  will  be  cooler  (denser)  than  its  surroundings  and 
subjected  to  a  downward  force.  Downward  moving  air  will  be  warmer  with 
an  upward  restoring  force.  In  either  case,  work  is  required  to  move  a 
volume  of  air  from  its  equilibrium  altitude,  so  there  is  a  positive 
change  in  the  potential  energy  of  the  displaced  air. 


If  we  assume  that  the  displacements  are  small  enough  so  that  the 
higher  order  effects  can  be  ignored  (i.e.,  that  the  lapse  rate  is  con¬ 
stant  over  the  displacement  distance),  then  the  restoring  force  is 
proportional  to  the  temperature  difference  between  the  displaced  parcel 
and  its  environment,  which  in  turn  is  proportional  to  the  displacement 
distance  (Az)  and  the  difference  between  the  ambient  lapse  rate  (Y)  and 
the  adiabatic  (Y  ).  The  corresponding  acceleration  (force  per  unit 
mass)  can  be  written  as  follows  (e.g.,  Holmboe  et  al  ,  19^8). 
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The  change  in  potential  energy  is  obtained  by  integrating  the  above 
expression  with  respect  to  z  over  the  complete  displacement  distance. 
This  results  in  the  following  expression  for  the  acquired  potential 
energy  per  unit  mass  (E): 
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where  g  is  the  acceleration  due  to  gravity  and  T  is  the  average  tempera¬ 
ture  (°K)  through  the  layer  encompassed  by  Az  . 


Steady-state  Lagrangian  plume  models  have  invoked  energy  considera¬ 
tions  in  the  concept  of  a  critical  dividing  streamline  height  in  order 
to  develop  simple  descriptions  of  plume  behavior  around  terrain  obsta¬ 
cles  in  a  stably  stratified  atmosphere  (e.g.,  Lavery  et  al  ,  1982). 
McNider  et  al  (1984)  attribute  this  concept  to  Hunt  and  Snyder  (1980). 
The  critical  streamline  height  (Hc)  is  the  height  at  which  the  kinetic 
energy  of  the  flew  exactly  matches  the  potential  energy  gained  by  rais¬ 
ing  the  air  to  the  top  of  the  obstacle.  Air  above  this  height  goes  over 
the  obstacle.  Below  the  critical  streamline  height,  the  air  goes  around 
the  obstacle,  providing  that  wind  speed  increases  with  height.  The 
governing  relationship  can  be  written 
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The  left  side  is  the  kinetic  energy  (per  unit  mass)  of  the  flow  at  the 
height  of  the  critical  streamline.  The  right  side  is  the  potential 
energy  (per  unit  mass)  gained  in  lifting  the  air  from  the  critical 
streamline  height  to  the  height  of  the  top  of  the  obstacle  (Z).  If  the 
vertical  temperature  gradient  is  constant  between  Hc  and  Z,  then  HQ  can 
be  expressed  in  terms  of  the  bulk  Froude  number: 


Hc  -  Z(1  -  Fr) 


where 
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For  our  purposes,  it  is  more  useful  to  restate  this  equivalence  of 
potential  and  kinetic  energy  in  terms  of  a  critical  height,  Zmax,  the 
highest  height  to  which  the  air  at  height  z  could  be  propelled  by  the 
local  wind  speed  (u)  against  the  local  30/3z.  That  critical  height  is 
given  by 


(6) 


III  MASS  CONSERVATION 


As  noted  earlier,  mass-conserving  wind  interpolation  schemes  obtain 
an  "initial  guess"  and  then  adjust  that  guess  to  remove  the  divergence. 
An  initial  guess  for  the  winds  above  the  surface  are  more  difficult  to 


An  initial  guess  for  the  winds  above  the  surface  are  more  difficult  to 
obtain.  Sherman  (1978)  used  the  shape  of  any  measured  wind  profile  (when 
available)  to  extrapolate  upward  throughout  the  region.  When  no 
observed  profile  was  available,  Sherman  assumed  a  linear  variation  of 
wind  with  height.  Endlich  et  al.  (1982)  assumed  a  uniform  wind  at  the 
top  of  the  domain.  They  derived  that  wind  from  the  surface  pressure 
field,  using  the  geostrophic  wind  relationship  (see  e.g.,  Holmboe, 
1 948 ) .  More  recent  versions  of  the  model  (Endlich  and  Lee,  1983)  use 
both  surface  temperature  and  pressure  observations  with  both  the  geo¬ 
strophic  and  thermal  wind  relationships  to  provide  upper  level  wind 
estimates.  The  most  recent  version  of  this  work  is  designed  to  use 
Doppler  acoustic  sounder  wind  profiles  with  the  geostrophic/thermal  wind 
estimate  to  obtain  a  spatially  varying  wind  field  at  the  top  of  the 
domain  (Nitz  et  al . ,  1985).  Winds  at  intermediate  levels  are  derived 
from  a  log-linear  interpolation  between  the  surface  wind  estimate  and 
that  for  the  top  of  the  domain. 

In  addition  to  the  differences  in  the  way  in  which  the  initial  wind 
field  estimate  is  obtained,  there  are  also  major  differences  in  the 
coordinate  systems  that  are  used.  For  example,  Sherman  (1978)  uses 
standard  rectangular  coordinates  while  Endlich  et  al.  (1982)  use  a 
curvilinear  coordinate  system  that  depends  on  the  shape  of  the  upper  and 
lower  boundaries  of  the  domain.  Ludwig  (1985)  has  reviewed  the  main 
features  of  these  two  modeling  approaches  and  discussed  some  of  the 
difficulties  of  incorporating  energy  conservation  schemes  into  the 
rectangular  coordinate  system.  The  following  discussion  is  limited  to 
curvilinear  coordinates. 

Bhumralkar  et  al.  (1980)  used  a  curvilinear  coordinate  system  to 
simplify  the  treatment  of  the  lower  boundary  in  their  mass-consistent 
wind  interpolation  scheme.  That  work  was  modified  by  Endlich  et  al. 
(1982)  and  later  by  Nitz  et  al .  (1985).  This  latter  version  (Nitz  et 
al.,  1985)  forms  the  basis  for  the  mass  conservation  aspects  of  the 
model  described  here.  In  this  model,  the  flow  is  assumed  to  take  place 
within  specified  surfaces.  We  will  refer  to  these  as  "flow"  surfaces. 
They  are  similar  in  many  respects  to  the  "sigma"  surfaces  used  by  others 
(e.g.  Bhumralkar  et  al,  1980),  but  they  differ  in  at  least  one  important 
respect.  The  flow  surfaces  can  intersect  the  underlying  terrain. 
Wherever  they  intersect  the  terrain,  the  wind  will  be  zero,  and  mass 
conservation  will  require  that  there  be  flow  around  the  obstacle.  The 
computational  scheme  allows  for  these  internal  boundaries. 

Before  proceeding  to  discuss  how  the  topography  of  the  flow  sur¬ 
faces  is  determined,  we  will  review  the  methods  by  which  non-divergence 
is  achieved  and  vertical  velocity  is  calculated.  Inasmuch  as  the  sur¬ 
faces  are  assumed  to  conform  with  the  streamlines,  vertical  velocity  (in 
a  rectangular  system)  is  given  by 
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where  is  the  horizontal  velocity  and  Vzsfc  is  the  slope  of  the  flow 
surface.  If  the  slope  is  not  too  great,  then  the  velocity  in  the  flow 
surface  can  be  used. 

Nitz  et  al  (1982)  replaced  the  wind  component  variables  with  the 
following 

u*  =  uAz 

(8) 

v*  =  vAz 

where  Az  is  the  average  separation  between  the  surface  and  the  next 
higher  and  next  lower  ones.  At  the  bottom  and  top  of  the  domain,  Az  is 
half  the  distance  to  the  adjacent  surface.  If  the  vertical  motions 
caused  by  the  slope  of  the  flow  surfaces  are  small  compared  to  the 
horizontal  motions,  the  continuity  equation  can  be  written 
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The  iterative  scheme  for  producing  nondivergence  has  been  described  by 
Endlich  (1984)  and  Nitz  et  al  (1985). 


IV  DEFINING  FLOW  SURFACES 

In  essence,  the  objective  is  to  define  a  coordinate  system  that  is 
more-or-less  flow  following.  Once  this  is  done,  then  there  is  no  reason 
why  the  coordinate  surfaces  may  not  intersect  the  terrain  in  areas  where 
there  is  flow  around  an  obstacle.  The  winds  are  set  to  zero  for  any 
grid  points  that  lie  below  the  local  terrain  height.  As  noted  before, 
this  approach  has  been  implemented  in  a  recent  version  of  the  model 
(Nitz  et  al.,  1985).  In  order  to  make  the  method  work  properly,  there 
must  be  some  way  to  determine  where  a  given  sigma  surface  intersects  the 
terrain.  Nitz  et  al.  (1985)  provide  only  subjective  guidance.  The 
critical  streamline  concept,  as  stated  in  Equation  (6)  seems  to  provide 
a  basis  for  an  objective  approach  to  the  determination  of  sigma  surface 
shapes.  That  equation  can  be  rewritten 


where  the  subscript  "o"  refers  to  conditions  at  the  lowest  geometric 
altitudes  found  on  a  particular  flow  surface.  These  conditions  are  then 
used  in  Equation  (10)  to  determine  zmax,  the  highest  altitude  that  that 
particular  flow  surface  will  reach.  In  order  to  implement  a  system 
based  on  this  approach,  the  following  steps  are  used: 
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(1)  The  lowest  and  highest  surface  elevations  within  the  model 
domain  are  identified. 

(2)  The  top  of  the  domain  (i.e.  the  top  of  the  boundary  layer) 
over  the  lowest  terrain  elevation  is  defined  and  the  vertical 
extent  of  the  domain  at  that  point  is  divided  into  uniformly 
separated  layers. 

(3)  A  value  of  V  is  estimated  from  the  initially  interpolated 
wind  field  for  each  flow  level;  the  estimate  is  derived  from 
an  average  of  the  interpolated  winds  at  the  appropriate  levels 
above  each  of  the  five  lowest  grid  points.  Other  parameters 
in  Equation  (10)  are  estimated  from  a  temperature  soinding, 
although  Doppler  acoustic  sounder  observations  of  ow  or  clima¬ 
tological  information  could  probably  be  used. 

(4)  Equation  (10)  is  solved  for  for  each  flow  surface;  these 

define  the  heights  of  the  flow  surfaces  at  the  x,y  coordinates 
of  the  highest  topographical  point. 

(5)  The  value  calculated  for  each  flow  surface  is  adjusted 

(beginning  at  the  second  from  topmost  level)  so  that  the  lower 
surfaces  do  not  approach  too  close  to  the  surfaces  above  them; 
surfaces  are  currently,  restricted  from  being  closer  than  half 
their  separation  above  the  lowest  terrain  elevations. 

(6)  A  second,  upward  pass  is  made  through  the  z[Dax  list  to  insure 
that  each  upper  surface  is  parallel  to  the  lower  surfaces,  or 
approaches  closer  to  them  over  the  high  terrain;  this  prevents 
the  effects  of  the  underlying  terrain  from  penetrating  upward 
through  stable  layers. 

(7)  The  heights  of  flow  surfaces  above  the  remaining  grid  points 
are  interpolated  linearly  according  to  the  surface  terrain 
he  i ght . 

(8)  Mass-consistent  wind  field  adjustments  are  made  using  a  model 
like  that  of  Nitz  et  al .  (1985). 

The  procedure  outlined  above  assumes  a  positive  (stable)  vertical 
potential  temperature  gradient  throughout  the  depth  of  the  modeling 
domain.  If  this  is  not  the  case,  then  the  flow  surfaces  will  be  terrain 
following,  i.e.  ,z-max  is  equal  to  the  difference  in  elevation  between 
highest  and  lowest  surface  points.  A  standard  interpolation  scheme  in 
rectangular  coordinates  is  applied  to  obtain  a  "first-guess"  wind  field 
that  is  used  for  the  flow  surface  definition  procedure  just  described. 
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FIGURE  C-T . 


Terrain  used  for  wind  model  test. 
(Contour  interval  is  5  meters) 


V  TEST  RUNS 


The  above  procedure  has  been  implemented  and  preliminary  tests  have 
been  completed.  A  21  x  21  x  6  grid,  extending  from  the  surface  to 
about  200  m  was  superimposed  on  the  terrain  shown  in  Figure  C-1.  Two 
cases  were  studied.  In  both  cases,  a  steady  south  wind  was  used  with 
wind  speeds  of  1  m  s"1  near  the  surface  and  2  m  s-1  at  an  altitude  of 
300  m  above  the  plain  on  which  the  92  m  hill  is  located.  This  wind 
profile  was  used  with  a  neutral,  adiabatic  lapse  rate  and  with  a  stable, 
inversion  rate  of  3. 6  x  10_2°C  m” 1 . 


Figure  C-2  shows  the  shape  of  lowest  flow  surface  for  the  two 
tests.  The  upper  part  of  the  figure  shows  that  the  flow  surface  paral¬ 
lels  the  topography  throughout  the  region  for  the  neutral  case.  The 
stable  stratification  has  a  modest  rise  as  the  terrain  is  approached, 
but  the  surface  does  not  pass  over  it.  This  forces  a  flow  around  the 
hill  as  can  be  seen  in  Figure  03.  For  the  neutral  case,  the  vectors 
are  all  parallel,  indicating  flow  over  the  hill.  Some  slowing  of  the 
wind  arises  over  the  elevations  near  ^0  m  because  of  the  logarithmic 
interpolation  upward  from  the  zero  wind  speed  at  the  surface.  The  flow 
around  the  hill  is  clearly  evident  in  the  stable  case  at  this  altitude. 


Figure  C— -4  shows  the  results  at  the  level  80  m  above  the  terrain 
surrounding  the  hill.  Again,  the  neutral  case  shows  wind  directions  to 
be  unaltered  by  the  hill.  The  speed  reductions  are  less  pronounced  than 
at  the  lower  altitude.  The  next  highest  flow  surface  (at  80  m  above  the 
flat  terrain)  just  barely  clears  the  top  of  the  hill  for  the  stable  case 
and  causes  deviations  in  the  flow  so  that  it  is  nearly  parallel  to  the 
topographical  contours  in  some  places.  There  is  even  some  flow  through 
the  small  col  at  the  crest  of  the  hill. 
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Figure  C-5  shows  two  cross  sections  through  the  center  of  the  hill 
from  south  to  north.  The  vertical  scale  is  magnified  by  a  factor  of  3 
relative  to  the  horizontal,  the  same  difference  in  scale  is  used  for  the 
w  (vertical)  and  v  (northward)  components  of  the  plotted  vectors.  The 
terrain  following  streamlines  are  clearly  evident  in  the  neutral  test 
results.  However,  the  stable  test  shows  streamlines  rising  only 
slightly. 
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FIGURE  C-5.  Flow  in  cross-sections  through  the  center  of  the  hill  for 
neutral  and  stable  test  cases. 


The  results  obtained  here  required  only  about  15  seconds  cpu  time 
on  a  VAX  11/782  computer.  The  storage  requirements  are  mostly  accounted 
for  by  7  floating  point  arrays  that  must  be  at  least  as  large  as  the 
grid  used  (21  x  21  x  6  in  this  case),  and  two-dimensional  arrays  (21  x 
21  for  these  tests).  There  are  also  numerous  smaller  arrays,  but  over¬ 
all  the  computational  requirements  are  limited  enough  that  the  technique 
should  be  usable  in  real  time  on  machines  that  are  small  enough  for 
routine  field  use. 

The  model  discussed  has  not  been  fully  tested.  We  hope  to  compare 
model  outputs  with  observations  from  regions  of  varying  topographical 


complexity,  and  to  use  the  results  to  refine  the  model.  The  model 
already  shows  considerable  promise  for  practical  applications.  The 
early  results  are  certainly  qualitatively  realistic.  The  data  require¬ 
ments  are  modest,  and  it  is  likely  that  a  single  Doppler  acoustic 
sounder  could  provide  the  information  necessary  to  run  the  model, 
although  more  inputs  would  improve  the  results. 

The  FORTRAN77  source  code  for  this  model  follows  the  references. 
Tables  C-1  through  C-3  summarize  the  input  data  file  requirements.  They 
also  show  the  logical  unit  designations  from  which  the  various  inputs 
are  read.  Outputs  (and  the  corresponding  logical  units  for  the  output 
files)  are  summarized  in  Tables  C-^  and  C-5.  Grid  sizes  are  defined  in 
data  statements  in  the  code. 

Most  of  the  input  variable  requirements  are  self-explanatory. 
However,  if  the  critical  streamline  feature  is  not  invoked,  the  user 
should  be  familiar  with  the  parameters  that  are  used  to  set  the  shape  of 
the  flow  surfaces. 

The  thickness  of  the  boundary  layer  and  the  shape  of  the  top 
(whether  curved  or  flat)  are  important  factors  in  producing  the  charac¬ 
teristics  of  the  nondivergent  flow.  In  many  instances,  radiosonde  data 
that  describe  the  boundary  layer  top  (BLT)  are  not  available.  In  such 
cases,  average  thickness  over  the  grid  (taken  from  experience  or  clima¬ 
tological  values)  and  a  slope  factor  are  used  to  specify  the  height  of 
the  BLT  at  each  grid  point.  The  height  of  the  BLT  at  a  particular  grid 
point  is  denoted  by  H;  the  average  thickness  of  the  boundary  layer  in 
the  area  of  interest  is  denoted  by  H;  hg  is  the  terrain  height  at  the 
grid  point;  hQ  is  the  terrain  height  at  a  reference  point  where  H  =  H  + 
hQ;  and  k  is  the  slope  factor.  Then  the  basic  equation  is 

H  =  H  +  khg  +  (1  -  k)hQ 

If  k  =  0,  the  boundary  layer  top  is  flat;  if  k  «  1,  the  top  parallels 
the  terrain.  Values  between  zero  and  one  give  intermediate  slopes,  and 
values  greater  than  one  give  slopes  steeper  than  the  terrain  slope. 
Negative  values  give  slopes  opposite  to  the  terrain  slopes.  Figure  C-6 
shows  typical  boundary  layer  tops  for  a  nighttime  case  (H  =  500  m,  k  = 
0.2)  and  a  daytime  case  (H  =  1500  m,  k  =  0.8).  The  parameters  H  and  k 
can  be  treated  as  functions  of  time  of  day  and  season. 

In  the  vicinity  of  coastlines  the  BLT  usually  has  a  general  slope 
from  ocean  to  higher  values  over  land  in  the  daytime,  with  the  reverse 
situation  at  night.  To  account  for  this,  two  additional  parameters  are 
used  to  specify  a  tilt  to  the  BLT.  One  parameter  specifies  the  west-to- 
east  difference  in  the  thickness  of  the  boundary  layer  across  the  grid, 
and  the  other  parameter  specifies  the  south-to-north  difference.  There 
is  also  a  parameter  for  the  minimum  thickness  of  the  boundary  layer  over 
isolated  mountain  peaks.  These  five  parameters  are  used  in  an  iterative 
solution  for  the  height  of  the  BLT  (Endlich  and  Lee,  1983). 


FIGURE  06.  Typical  configurations  of  the  boundary- layer  top  for 
daytime  and  nighttime. 


Table  C-1 


TOPOGRAPHIC  ELEVATION  INPUTS 
(Logical  Unit  1 1) 


Format  _ Variables _ 

8F10.2  ( ( HT CRS ( I ,  JR ) ,  I-1.MCRS),  JR=NCRS,  1 1 ) 

8F10.2  ((HTMED(I.JR),  I-1.MMED),  JR-NMED , 1 1 ) 

8F10.2  ( (HTFIN(  I,  JR),  I=1,MFIN),  JR-NFIN, 1 ,-1 ) 


HTCRS,  HT  ME  D,HTF  IN --elevations  of  points  in  coarse, 
medium  and  fine  grids.  HTMED  read  only  if  more 
than  one  grid  used.  HTFIN  read  only  if  more  than 
two  grids  used. 

MCRS,  NCRS — number  of  coarse  grid  points  in  x,y 
directions  (=21  x  21). 

MMED.NMED — number  of  medium  grid  points  (21  x  21 )  . 
MFIN.NFIN — number  of  fine  grid  points  (21  x  21). 


Table  C-2 


VARIABLES  READ  BY  WIND  MODEL  FROM  LOGICAL  UNIT  12 


Format 

Variabl  es 

Subroutine 

* 

AUTHK,  SLFAC,  STNK ,  DNI,  BLGRX,  BLGRY 

MAIN 

# 

TERLIM,  IDATE,  IHOUR 

MAIN 

# 

( LPRNT(K) ,  K-1,  NLVL 

MAIN 

* 

NUMDOP 

DOPSIG 

# 

(STAID(IT),  STLT(IT),  STLN(IT),  IT=1, NUMDOP) 

DOPSI G 

* 

(NHTS,  continued 

* 

(DPHT( IT, LL) ,  DPWD ( IT , LL) ,  DPWS(IT,LL),  LL=1 , NHTS) , 

DOPSIG 

ISITE-1 , NUMDOP) 

* 

NUMNWS 

WXANAL 

# 

(STAID(IT),  STLT(IT),  STLN( IT) ,  PRESS(IT),  TEMP( IT) , 

GEOSI G 

WD(IT),  WS(IT),  IT-1, NUMNWS) 

# 

NTSOND 

WXANAL 

AVTHK  =  average  boundary  layer  thickness  (m) 

SLFAC  =  degree  to  which  terrain  is  followed  0=flat,  1=follow 
BLGRX,  BLGRY  =  boundary  layer  gradient 
TERLIM  -  intersection  height  with  terrain 

The  above  parameters  used  primarily  when  the  stability  option 
is  not  used.  AVTHK  defines  heights  of  surfaces  over  low  terrain, 
otherwise  the  values  are  overridden  by  stability  method. 

IDATE,  IHOUR  «  date  and  hour  for  identification 
LPRNT  =  levels  to  be  printed  have  LPRNT-1 
NUMDOP  =  number  of  Doppler  sodars 
STAID  =  Doppler  station  ID  numbers  (not  used) 

STLT,  STLN  =  Doppler  station  latitude  and  longitude  (decimal  degree) 
NHTS  =  number  of  elevations  for  Doppler  input 

DPHT ,  DPWD,  DP WS  =  height  (m),  wind  direction  (deg),  speed  (m  s-1) 
NUMWS  =  number  of  sfc  wx  inputs 

STAID  =  sfc  station  identification  number  (not  used) 

STLT,  STLN  =  sfc  station  latitude  and  longitude 

NTS ON D  =  indicates  if  stability  option  used  (=1)  or  not  used  (=0) 

Note  Units: 

PRESS,  TEMP  =  sfc  station  pressure  (”Hg)  and  temperature  (° F) 

WD.WS  =  sfc  station  wind  direction  (decimal  deg)  and  speed  (m  s-1) 
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Table  C-3 


TEMPERATURE  SOUNDING  INPUTS 
(LOGICAL  UNIT  13) 


Format  _ _ Variables* _ 

*  NHITES ,  ITYP 

*  (2(1),  T(I),  P(I),  1=1 , NHIES) 


NHITES  —  number  elevations  for  inputs 
ITYP  —  type  of  inputs  (see  below) 

Z,T,P  —  height  (m),  temperature  (°C), 
de/dz  for  ITYP-1 

Z,T,P  —  height  (m),  temperature  (°C), 
pressure  (mb)  for  ITYP=0 


*Read  only  if  NTSOND  (Table  C-2)  is  1. 


Subroutine 

STRAT 

STRAT 


Table  C-4 

WIND  COMPONENT  OUTPUTS 
(LOGICAL  UNIT  33) 

Format 

Variables 

IX, 2116  ( (U(I,J,L) ,  1-1,  NCOL) ,  J-NROW, 1 , 1-1 ) 

IX, 2116  ( (V( I , J,L ) ,  1-1,  NCOL),  J-NROW, 1,-1) 

IX, 2116  ( (W( I , J , L) ,  1=1,  NCOL),  J-NROW, 1,-1) 

REPEATED  FOR  L-2.NLVL 


U,V,W  -  wind  components  (cm  s~' ) 


NROW,  NCOL,  NLVL  =  number  of  rows,  columns,  and 
levels 


Winds  are  for  flat  surfaces,  at  heights  determined 
by  flow  surface  heights  above  the  lowest  terrain. 


Table  C-5 


FLOW  SURFACE  TOPOGRAPHY  OUTPUTS 
(LOGICAL  UNIT  15) 


Format 


(1X.25F5.0) 


Variables 


(((RHS(IX,IY,L)  +  SFCHT ( IX , IY ) ,  IX-1,21) 
IY-21  ,1 ,-1 ) ,  L=1 ,  6) 


RHS  =  height  of  flow  surface  relative  to  underlying 
terrain  (m) 


3FCHT  =  height  of  underlying  terrain  (m) 


•  ’  ■  *  »  *  »  *  «"*. 
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PROGRAM 


WND 


••THIS  FORM  OF  PROGRAM  COMPLEX  COMPUTES  TOPOGRAPHICALLY  INDUCED 

*  WINDS  BY  MAKING  THE  ORIGINALLY  ANALYZED  WINDS  NONDIVERGENT . 

SIGMA  SURFACES  ARE  DEFINED  USING  A  CONCEPT  ANALAGOUS  TO 
THE  CRITICAL  STREAMLINE  WHEN  THE  ATMOS  IS  STABLY  STRATIFIED. 

THESE  MODIFICATIONS  ARE  PRIMARILY  CONTAINED  IN  SUBROUTINES 
XTREM, STRAT,k£SIG .  SIGMA  SURFACES  CAN  INTERSECT  THE  TERRAIN 
RESULTING  IN  ZERO  WINDS  IN  THE  "UNDERGROUND  CELLS" 

THIS  FORM  ALSO  REPLACES  INWND  WITH  WXANAL  WHICH  READS  IN  WIND 
SOUNDINGS  AND  NWS  DATA  AND  MAKES  THE  INITIAL  WIND  ANALYSIS. 

*  THIS  FORM  USES  DIRECT  VECTOR  ALTERATIONS  IN  SUBROUTINE  BAL5 . 

OTHER  NEW  FEATURES  INTRODUCED  IN  '84-' 85  ARE: 

WIND  COMPONENTS  AT  SPECIFIED  POINTS  CAN  BE  HELD  CONSTANT  IN  BAL5 . 
BY  R.  ENDLICH,  F.  LUDWIG,  K.  NITZ,  C.  MAXWELL 
SRI  INTERNATIONAL 

333  RAVENSWOOD  AVE  MENLO  PARK,  CALIF  94025. 

PHONE  (415)  859-3395  OR  (415)  859-2915 
INTEGER  DNI 

LOGICAL  DBUG (15),  KEY,  UTMUSE 

DIMENSION  B( 25 , 25 ) , LPRNT( 10 ) , IB( 25 , 25 ) , ZFLAT( 6 ) 

COMMON/RARS/RHS ( 25 , 25 , 6 ) 

COMMON/CSFC/SFCHT (25,25), SIGMA ( 6 ) 

COMMON /UARS/U( 25,25,6),UA(25,25,6),V(25,25,6),VA(25,25,6) 

COMMON /WARS/W( 25 , 25 , 6 ) , WA( 25 , 25 , 6 ) 

COMMON/PARMS/ZTOP , DS , DS IGMA , NLVLM1 , XHT1 , XHT2 , XI , Y1 , 

1  X2,Y2,UG,VG, RATIO, TDSI 

COMMON  /CVOS/  RCM,RMF, IV,DSCRS, IXCRS, JYCRS, IXMED, JYMED, 

+  IXFIN , JYFIN 

COMMON/CTOP/  MCRS , NCRS , MMED , NMED , MFIN , NFIN , NGR ID 
COMMON  /LIMITS/NCOL, NROW, NLVL, NCOLM1 , NROWM1 , 

$  LOWIX( 5,3), LOWIY (5,3), SFCMAX 

COMMON  /BLHT/  BLT( 25 , 25 ) , HSITE,  AVTHK ,  SLFAC , STHK , BLGRX , BLGRY 

COMMON  /SITE/  IXSA  JYS,  THSITE,  IGRID 

COMMON  /ANCHOR/  SLAT,  SLNG  , UTMUSE 

COMMON  /GROUND/  TERLIM 

COMMON  /PRL1M/,  II,  12 

C  SET  LIMIT  ON  NO.  OF  ITERATIONS  IN  SUBR .  BAL5 
DATA  NIT/15/ 

C  SET  NUMBER  OF  VERTICAL  LEVELS 
DATA  NLVL/6/ 

C  SET  SIGMA  LEVELS  (VERTICAL  COORDINATES) 

DATA  SIGMA/0.0,  0.2,  0.4,  0.6,  0.8,  1.0/ 

C  SET  NO.  OF  GRIDS  TO  BE  USED 
DATA  NGR ID/1/ 

C  SET  NO.  COLS, ROWS  FOR  COARSE , MED , FINE  GRIDS 

DATA  MCRS/ 17/ , NCRS/ 13/ , MMED/ 21/ , NMED/ 21/ 

DATA  MFIN/21/ , NFIN/21/ 

C  SET  GRID  INTERVALS  IN  KM 

DATA  DSCRS/0 . 5/ ,  DSMED/0.25/,  DSFIN/0 . 1/ 

C  SET  RATIOS  OF  GRID  INTERVALS,  COARSE  TO  MED.  ,  MED.  TO  FINE 
DATA  RCM/2 . 0/ ,  RMF/2 . 5/ 

C  SET  HEIGHT  (METERS)  OF  ANCHOR  POINT  (REFERENCE  POINT) 

DATA  HSITE/838 . 0/ 

C  SET  LATITUDE,  LONGITUDE  OF  ANCHOR  POINT  ( UTMUSE= . FALSE ) ,  OR 
C  UTM  COORDINATES  OF  ANCHOR  PT  ( UTMUSE= . TRUE . ) 

DATA  SLAT/4290.0/,  SLNG/522.0/  , UTMUSE/ . TRUE . / 

C  SET  ANCHOR  POINT  LOCATION  (IX,JY)  IN  EACH  GRID 

DATA  IXCRS/  1/, JYCRS/  1/ , IXMED/ 7/ , JYMED/  5/ 

DATA  IXFIN/ 7/ , JYFIN/ 5/ 

C  SET  LIMITS  FOR  PRINTING  COLUMNS 

DATA  I 1/1/,  12/17/  , KEY/. TRUE./ 

C 

C  DBUG  CALLS  PRINTS  IN  BAL5  DOPSIG  GEOSIG  GPAN  RESIG 
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n a  non  no  noon  no  on 


DATA  DBUG  /.FALSE.,  .FALSE.,  .FALSE.,  .FALSE.,  .FALSE., 

SETBLT  SETMAT  STRAT  TOPO  WXANAL 

$  .FALSE.,  .FALSE.,  .FALSE.,  .FALSE.,  .FALSE., 

XTREM  MAIN  PROG  REINT  LEVWND  1ST  ESTIMATES 
$  .FALSE.,  .TRUE.,  .FALSE.,  .FALSE.,  .TRUE./ 

*  IN  ALL  ARRAYS  POINT  (1,1)  IS  AT  SW  CORNER.  X  INCREASES  TO  EAST,  Y 

TO  N.  INDICES  ARE  I,J,K  -(COL, ROW, LYR )  WITH  LIMITS  NCOL , NROW , NLVL 
UNITS  USED  IN  COMPUTATIONS  ARE  M,  G,  SECONDS 

IF(DBUG( 12 ) )  PRINT  9013,  MCRS, NCRS , MMED, NMED 
IF ( DBUG ( 12 ) )  PRINT  9014,  DSCRS,  DSMED,  DSFIN 
READ  BOUNDARY  LAYER  VARIABLES  (USED  IN  SETBLT) 

READ( 12 , 9022 )  AVTHK ,  SLFAC ,  STHK,  DNI,  BLGRX ,  BLGRY 
READ (12,*)  AVTHK,  SLFAC,  STHK,  DNI,  BLGRX,  BLGRY 
IF ( DBUG ( 12 ) )  PRINT  9025,  AVTHK,  SLFAC,  STHK,  DNI,  BLGRX,  BLGRY 
UNITS  ARE  METERS 

READ  TERRAIN  INTERSECTION  (USED  IN  TOPO),  DATE,  HOUR 
READ  (12,9023)  TERLIM,  IDATE,  IHOUR 
READ  (12,*)  TERLIM,  IDATE,  IHOUR 
IF(DBUG( 12 ) )  PRINT  9026,  TERLIM 
IF(DBUG( 12 ) )  PRINT  9027,  IDATE,  IHOUR 

*  READ  PRINT  INDICATORS.  TO  PRINT  LEVEL  K  USE  LPRNT(K)=1. 

READ( 12 , 9030 )  (LPRNT(K) ,K=1,NLVL) 

READ( 12 , * )  (LPRNT(K) ,K=1,NLVL) 

IF ( DBUG ( 12 ) )  PRINT  9035 

I F ( DBUG (12))  PRINT  9030,  ( LPRNT( K ) , K=1 , NLVL ) 

C  LOOP  THRU  NGRID  SYSTEMS. 

DO  1040  IGRID=1, NGRID 
IF( DBUG ( 12 ) )  THEN 

IF  ( IGRID  .EQ.  2)  PRINT  9005 
IF  (IGRID  .EQ.  3)  PRINT  9006 

END  IF 

C  SET  CONTROLS  FOR  PROPER  GRID 

IF  (IGRID  .NE.  1)  GO  TO  11 
IXS=IXCRS 
JYS=JYCRS 
NCOL=MCRS 
NROW=NCRS 
DS=DSCRS 
11  CONTINUE 

IF  (IGRID  .NE.  2)  GO  TO  14 

IXS= IXMED 

JYS= JYMED 

NCOL=MMED 

NROW=NMED 

DS= DSMED 

14  CONTINUE 

IF  (IGRID  .NE.  3)  GO  TO  15 

IXS=IXFIN 

JYS=JYFIN 

NCOL=MFIN 

NROW=NFIN 

DS=DSFIN 

15  CONTINUE 

C  READ  TERRAIN  HEIGHTS  FOR  ALL  GRIDS  USING  TOPO 
NUM1=0 

IF  (IGRID  .EQ.  1)  CALL  TOPO ( NUM1 , DBUG ) 

DS=DS*1 . 0E3 
NCOLMl=NCOL-l 
NROWM1 =NROW- 1 
NLVLM1=NLVL-1 
TDSI  =  1 ./( 2 . 0  *DS ) 

<:  SET  BL  TOP  AND  SIGMA  SURFACES 
CALL  TOPO ( IGRID, DBUG) 


C  PRINT  +  PLOT  SURFACE  HEIGHT 
I F ( DBUG (15))  PRINT  171 
DO  53  JP=1 , NROW 
DO  53  IP=1 , NCOL 

B{ IP, JP)=  SFCHT{ IP, JP) 

IB( IP, JP)=JNINT(B{ IP, JP) ) 

53  CONTINUE 

IF  ( DBUG (15))  THEN 
DO  54  JP=1 , NROW 

JR=NROW+  1-JP 

54  PRINT  9105,  ( IB( IP , JR ) , IP= II , 12 ) 

END  IF 

C  *  *  *  *PLOT  GEOMETRIC  HEIGHTS  OF  SELECTED  SIGMA  SURFACES 
DO  511  K=1 , NLVL 

IF  (LPRNT(K)  .NE.  1)  GO  TO  511 
DO  627  JP=1 , NROW 
DO  627  IP=1 , NCOL 

B(  IP, JP)=  RHS( IP, JP, K) 

IB( IP , JP ) =JNINT( B{ IP,JP) ) 

627  CONTINUE 

I F ( DBUG (15))  PRINT  571, K 
DO  628  JP=1 , NROW 
JR=NROW+  1-JP 

628  PRINT  9105,  ( IB( IP , JR ) , IP= 1 1 , 1 2 ) 

511  CONTINUE 

C*  READ  AND  ANALYZE  WIND  DATA  USING  WXANAL 
C*  MAKE  INITIAL  WIND  ANALYSIS  ON  MESH 
CALL  WXANAL ( IGR ID, DBUG) 

C  »  *  *  *  *  ‘PLOT  OBSERVED  VELOCITY  COMPONENTS  AT  SELECTED  LEVELS 
DO  211  K  =  1 , NLVL 
IF  (LPRNT(K)  .NE.  1)  GO  TO  211 
DO  46  J  =  1 , NROW 
DO  46  1=1, NCOL 

IB( I , J ) = JN1NT( U ( I , J , K ) *  10 .  ) 

46  CONTINUE 

IF( DBUG ( 15 ) )  PRINT  271, K 
DO  56  JP=1 , NROW 
JR=NROW+  1-JP 

56  IF(  DBUG(  15  )  )  PJUNT9105,  ( IB(  I ,  JR  )  ,  1  =  11, 12) 

211  CONTINUE 

DO  212  K=1 , NLVL 

IF  (LPRNT(K)  .NE.  1 )  GO  TO  212 
DO  48  J  =  1 , NROW 
DO  48  1=1, NCOL 

IB( I , J ) =JNINT( V ( I , J , K ) *10 . ) 

48  CONTINUE 

I F ( DBUG (15))  PRINT  272, K 
DO  58  JP= 1 , NROW 
JR=NROWf  1-JP 

58  I F ( DBUG (15))  PRINT  9105,  ( 1B( I , JR ) , I  =  1 1 , 1 2 ) 

212  CONTINUE 

IF  (IGRID.GT.l)  GO  TO  226 
C  COMPUTE  VERTICAL  MOTION  W  ALONG  SIGMA  SURFACES 
DO  220  K=  2 , NLVL 
DO  220  1=2, NCOLM1 

DO  220  J  =  2 , NROWM1 
W(  I , J , K )  =  0 . 0 

IF  ( RHS( I,J,K) .LE.0.0)  GO  TO  220  !  FOR  TERRAIN  LIMIT 

HSIGE= SFCHT (I+1,J)+RHS(I+1,J,K) 

HSIGW=SFCHT ( I- 1 , J ) +RHS ( I - 1 , J , K ) 

HSIGN=SFCHT{ I , J+ 1 ) +RHS ( I , J+ 1 , K ) 

HSIGS=  SFCHT { I , J-l )+RHS( I, J-l, K) 

DHDX= ( HSIGE- HS IGW ) *TDSI 
DHDY= ( HSIGN-HSIGS ) *TDSI 
W(  I  , J, K)=U( I , J, K) *DHDX+V( I , J, K ) *DHDY 
42U  CONTINUE 


n  o  n 


2  26  CONTINUE 

IF  ( IGRID  .Eg.  1)  PRINT  9050 
IF  (IGRID  .EQ.  2)  PRINT  9055 
IF  (IGRID  .Eg.  3)  PRINT  9056 
I F ( DBUG (15))  PRINT  9020 

I F ( DBUG (15))  PRINT  3,(K,  U( IXS , JYS , K ) ,  V( IXS, JYS, K ) , 
$  W( IXS, JYS, K ) ,  RHS( IXS, JYS, K) ,K=1, NLVL) 

DO  213  K=1 , NLVL 

IF  (LPRNT(K)  .NE.  1 )  GO  TO  213 
I F ( DBUG (15))  PRINT  273, K 
DO  60  JP=1 , NROW 
DO  60  IP=1 , NCOL 
IB( IP, JP)=JNINT(W( IP, JP, K) *100. ) 

60  CONTINUE 

I F ( DBUG (15))  THEN 
DO  65  JP=1 , NROW 
JR=NROW+  1-JP 

65  PRINT  9105,  ( IB( IP, JR } , IP= I 1 , 12 ) 

END  IF 

213  CONTINUE 

CALL  SUBROUTINE  TO  MAKE  WINDS  NONDIVERGENT 


CALL  BAL5( NIT, DBUG) 

CALL  LEVWND(DBUG, IGRID, ZFLAT) 

C 

C  PLOT  ADJUSTED  VELOCITY  COMPONENTS  INTERPOLATED  TO  SELECTED  LEVELS 

C 

DO  611  K=1 , NLVL 

IF  (LPRNT(K)  .NE.  1)  GO  TO  611 
DO  615  J=1 , NROW 
DO  615  1=1, NCOL 

IB( I , J ) = JN INT ( UA( I , J , K ) *  1 0 .  ) 

615  CONTINUE 

I F ( DBUG ( 12) ) THEN 
PRINT  671, K 
DO  622  JP= 1 , NROW 

JR=NROW+  1-JP 

622  PRINT  9105,  ( IB( I , JR ) , I  =  1 1 , 1 2 ) 

END  IF 

ell  CONTINUE 

DO  612  K=1 , NLVL 

IF  (LPRNT(K)  .NE.  1 )  GO  TO  612 
DO  616  J=1 , NROW 
DO  616  1=1, NCOL 

IB( I , J ) =JNINT( VA( I , J , K ) *  10 .  ) 

616  CONTINUE 

IF ( DBUG ( 12 ) )THEN 
PRINT  672, K 
DO  624  JP=1 , NROW 
JR=NROW+  1-JP 

c24  PRINT  9105,  ( IB( I , JR ) , I  =  1 1 , 1 2 ) 

END  IF 

6  I  2  CONTINUE 

DO  641  K=1 , NLVL 

IF  (LPRNT(K)  .NE.  1)  GO  TO  641 
DO  645  J=1 , NROW 
DO  645  1=1, NCOL 

IB( I , J ) =JNINT( WA( I , J , K ) *  10 .  ) 

645  CONTINUE 

I F ( DBUG (12) )THEN 
PRINT  675, K 
DO  652  JP=1 , NROW 

JR=NROW+  1-JP 


652 


PRINT  9105,  ( IB ( I , JR ) , 1  =  1 1 , I  2 ) 


END  IF 


n  o  o 


641  CONTINUE 

IF  (IGRID  .EQ.  1)  PRINT  9050 
IF  (IGRID  .EQ.  2)  PRINT  9055 
IF  (IGRID  .EQ.  3)  PRINT  9056 
IF( DBUG ( 12 ) )  PRINT  2 

I F( DBUG ( 12 ) )  PRINT  10, (K,UA( IXS, JYS , K ) ,VA(IXS,JYS,K) , 

$  RHS ( IXS , JYS , K ) , K=1 , NLVL ) 

C-- WRITE  WNDS  INTERPOLATED  TO  FLAT  SFCS  (OR  10M  ABOVE  TERRAIN  FOR  1ST  LVL) 

C 

DO  700  K=1 ,  NLVL 
DO  690  J=l,  NROW 
JR=NROW+l-J 

690  WRITE  (3,9065)  ( N1NT( 100 . *UA( I , JR , K ) ) , 1=1 , NCOL) 

DO  695  J=l,  NROW 
JR=NROW+l-J 

695  WRITE  (3,9065)  ( NINT( 100 . 0‘VA( I , JR , K) ) , 1=1 , NCOL) 

DO  697  J=l,  NROW 
JR=NROW+l-J 

697  WRITE  (3,9065)  ( N INT( 100 . 0 *WA( I , JR , K ) ) , 1=1 , NCOL ) 

700  CONTINUE 

C 

C - WRITE  WINDS  ON  SIGMA  SFCS - 

C 

C  DO  750  K=  2 ,  NLVL 

DO  990  J=l,  NROW 
C  JR=NROW+l-J 

C990  WRITE  (3,9065)  ( NINT( 100 . *U( I , JR , K ) ) , 1=1 , NCOL) 

C  DO  995  J=l,  NROW 

C  JR=NROW+l-J 

995  WRITE  (3,9065)  ( NINT( 100 . 0*V ( I , JR , K ) ) , 1=1 , NCOL ) 

DO  997  J=l,  NROW 
JR=NROW+l-J 

C997  WRITE  (3,9065)  ( NINT( 100 . 0*W( I , JR , K ) ) , 1=1 , NCOL ) 

C750  CONTINUE 

C - - - 

1040  CONTINUE 

2  FORMAT  (/’  FINAL  RESULTS  AT  ANCHOR  PT’/2X,’K  UA  VA 

+  REL.  HT§ . 1 / ) 

3  FORMAT ( IX, 1 2 , 2X , 2F10 . 2 , Fll . 3 ,  F10 . 1 ) 

10  FORMAT  (1X,I2,2X,2F10.2,F10.1) 

171  FORMAT  ( 1  HI ,  1  TERRAIN  HEIGHT,  METERS,  NORTH  ROW  FIRST1/) 

271  FORMAT ( 1 H 1 , '  U  COMPONENT  DECIMETERS/SEC,  LVL  =  '14/) 

272  FORMAT ( 1H1 , '  V  COMPONENT  DECIMETERS/SEC,  LVL  =  '14/) 

273  FORMAT ( 1H1 , '  W,  CM/SEC  ,  LVL  = ' 13/ ) 

571  FORMAT  ( 1H1 , '  HEIGHT  ABOVE  TERRAIN,  M,  LVL='I3/) 

671  FORMAT  (1H1,'  ADJUSTED  U  COMPONENT,  ', 

$  ’DECIMETERS/SEC,  LVL='I3/) 

672  FORMAT  ( 1H1 , '  ADJUSTED  V  COMPONENT,’, 

$  '  DECIMETERS/SEC,  LVL='I3/) 

675  FORMAT  (1H1,'  ADJUSTED  W  COMPONENT,', 

$  '  DECIMETERS/SEC,  LVL=’I3/) 

9005  FORMAT  (1H1,'  **BEGIN  COMPUTATIONS  FOR  MEDIUM  GRID**'/) 

9006  FORMAT  ( 1H1 ,  ’  *  *  BEGIN  COMPUTATIONS  FOR  FINE  GRID**'/) 

9013  FORMAT  (/'  COARSE  GRID  E-W'15,'  S-N'15,'  MEDIUM  GRID  E-W'14,'  S 
♦  ’14/) 

9 U 1 4  FORMAT  (/'  GRID  INCREMENTS  IN  KM,  COARSE= ' F4 . 1 ,  ’  MED.  =  'F4.1,' 

+  FINE= ' F4 . 1/ ) 

9020  FORMAT  (/'  ORIGINAL  U,  V,  W,  REL.  HTS  AT  ANCHOR  PT . ' / ) 

9022  FORMAT  ( F10 . 1 , F10 . 2 , F10 . 1 , 1 5 , 2F7 . 1 ) 

9023  FORMAT  (F10. 1,218) 

9o2  5  FORMAT  (//'  AVER.  BNDY .  THICKNESS  IN  M  = ' F8 . 1 , 

+  'SLOPE  FACTOR  FOR  BL  TOP= ’ F4 . 1 , ’  MIN.  THICKNESS= 1  F7  . 1 , 

+  /’  DAY1  NITE2  INDICATOR= ’ 1 3 ,  ’  B  LYR  GRADIENT  TO  EAST,  M='F7.1 

+  ’  TO  NORTH  = 1  F7 . 1/) 

9026  FORMAT  (/'  TERRAIN  INTERSECTION  WILL  OCCUR  FOR  HTS 


+  ABOVE 1 F6 . 1 , '  METERS ' / ) 

9027  FORMAT  (/'  DATE  ='I8,'  HOUR  = ' 14/ ) 

9030  FORMAT  (1215) 

9035  FORMAT  (/'  INDICATORS,  LPRNT{ K ) , FOR  PRINTING  FIELDS'/) 

9046  FORMAT  (1H1,'  V  COMP.  AT  LEVEL  3  '/) 

9050  FORMAT  {1H1,'  COARSE  GRID') 

9055  FORMAT  (1H1,'  MEDIUM  GRID1) 

9056  FORMAT  ( 1H1 , '  FINE  GRID' ) 

9060  FORMAT  ( 15 , E10 . 1 , 15 ) 

9065  FORMAT  (IX, 2116) 

9100  FORMAT  (/5X,22F5.0) 

9105  FORMAT  (/IX, 2515) 

STOP 

END 

C 

r»***«**«**********»»»»*****************»****»*************************** 

C 

SUBROUTINE  BAL5 ( NITER , DBUG ) 

C*  ****THIS  IS  VERSION  11,  OCTOBER  '84. 

C  THIS  ROUTINE  BALANCES  DIVERGENCE  TO  VALUES  IN  ARRAY  DD  (OR 
C  TO  DDIJ=0)  AND  VORTICITY  TO  ARRAY  VT(I,J). 

C  FOR  WIND  SPEED  IN  MPS .  DIV  AND  VORT  ARE  SCALED  TO  UNITS 
C  10  -6  SEC-1.  THE  METHOD  USES  DIRECT  VECTOR  ALTERATIONS. 

C  THIS  FORM  FOR  RECTANG  GRID  OMITS  TRIG  FUNCTIONS. 

C  THE  FLUX  FORMULATION  IS  USED  FOR  FINITE  DIFFERENCES.  FOR  SIGMA  LAYERS 
C  COMPUTE  NON-DIV  WINDS  FOR  WINDS  WEIGHTED  BY  THICKNESS  OF  LAYER. 

C  ASSUME  SIGMADOT=0 .  INDICES  IN  ARRAYS  (I,J,K)  ARE  I=COLUMN, 

C  J=ROW,  K= LEVEL;  PT  (1,1,1)  IS  AT  SW  CORNER  AT  GROUND. 

C  FOR  COMPUTATION  BOXES,  INDICES  REFER  TO  SW  CORNER  OF  BOX. 

C  I VORT  CONTROLS  USE  OF  VORTICITY.  IF=2  VORT  IS  NOT  HELD 
C  CONSTANT . 

C  BY  R.M.  ENDLICH,  SRI  INTN'L,  1ST  VERSION  FOR  LAYERS  JUNE  '82. 

C 

LOGICAL  DBUG (15) ,RYT3 

DIMENSION  VT(25,25) ,DI( 25,25 ),VO( 25,25 ) ,U1( 25,25 ), VI (25, 25) 
DIMENSION  UN(25,25),VN(25,25), THK( 25,25), NEWLVL( 5 ) 

DIMENSION  IFXPT( 25,25) 

COMMON  /UARS/U( 25,25,6) , UA( 25,25,6),V(25,25,6), VA( 25,25,6) 
COMMON/PARMS/ZTOP , DS , DSIGMA , NLVLM1 , XHT1 , XHT2 , XI , Y1 , 

+  X2,Y2,UG,VG, RATIO, TDSI 

COMMON  /LIMITS/NCOL , NROW , NLVL , NCOLM1 , NROWMl , 

$  LOWIX(5,3) ,LOWIY(5,3) , SFCMAX 

COMMON  /SITE/  IXS,  JYS,  THSITE,  IGRID 
COMMON  /RARS/RHS (25,25,6) 

COMMON  /PRLIM/  11,12 
DATA  NEWLVL/2 , 6 , 3 , 4 , 5/ 

DATA  RYT3/. FALSE./ 

IVORT=2  !  IGNORE  VORTICITY 
IF(DBUG( 1 ) )  PRINT  9002 
GS=DS*1 . 0E-05 

C  USE  GRID  SPACING  IN  100'S  OF  KM.  DS  IS  IN  M. 

GSI=10 . 0/GS  !  FOR  PROPER  SCALING 

C  FOR  RECTANGULAR  GRID  OMIT  TRIG  FUNCTIONS  USED  PREVIOUSLY 
C  ASSIGN  PTS  WHERE  INITIAL  WIND  ANALYSIS  IS  HELD  FIXED 
IF  (IGRID  .NE.  1)  GO  TO  10  !  FOR  COARSE  GRID 

CALL  SETMAT( 0 ,  IFXPT,  NCOL,  NROW) 

C 

C - IDENTIFY  PTS  NOT  TO  BE  CHANGED - 

C  IFXPT( 11 , 8 ) =1 

10  CONTINUE 

IF  (IGRID  .NE.  2)  GO  TO  15  !  FOR  MEDIUM  GRID 

CALL  SETMAT ( 0 ,  IFXPT,  NCOL,  NROW) 

C  IFXPT(  ,  )=  1 

15  CONTINUE 


0  PRINT  POINTS  WITH  FIXED  WINDS 
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C  DO  25  J=l,NROW 

C  DO  25  1=1 , NCOL 

C  IF  ( IFXPT( I ,  J )  .EQ.  1)  PRINT  9582,  I,J 

C  25  CONTINUE 


C  ARRAYS  U , V, RHS  ARE  WRITTEN  { COLS , ROWS , LEVELS ) 

C  TRANSFER  WINDS  TO  2D  ARRAYS  IN  ORDER  BOTTOM( LEV2 ) ,  TOP, 3, 4, 5 

C 

C  DO  800  NORDR=l , NLVL-1 

C 

C  INTERPOLATING  IN  VERTICAL  TO  GET  1ST  GUESS  FIELDS  FOR  LEVELS  3,4,5 
C 

C  L=NEWLVL( NORDR ) 

C  IF  (NORDR/GE.3)  CALL  REINT( DBUG , L) 


C  NOT  REINTERPOLATING  FROM  PREVIOUSLY  BALANCED  LEVELS 
C 

DO  800  L=2 , NLVL 


DO  3  5  J  =  1 , NROW 
DO  35  1=1, NCOL 

UN ( I , J ) =U ( I , J  ,  L ) 

VN( I , J)=V( I , J,  L) 

35  CONTINUE 

CALL  SETMAT ( 0 . 0 , DI  , NCOL, NROW) 

CALL  SETMAT (0.0, VT  , NCOL, NROW) 

CALL  SETMAT (0.0, VO  , NCOL, NROW) 

C  COMPUTE  LAYER  THICKNESS  AND  MULTIPLY  WIND  COMPONENTS 
DO  40  J  -1 , NROW 
DO  40  1=1, NCOL 
LA=L  +1 

IF  (LA.GT.NLVL)  LA=NLVL 
HTA=RHS ( I , J , LA ) 

LB=L-1 

IF  (LB.LT.l)  LB=1 
HTB=RHS( I, J,LB) 

IF  (HTB.LT.-1.0)  HTB=~1.0 
THK( I ,  J )=0 . 5* ( HTA-HTB) *0 . 01 


IF  ( THK ( I , J ) . LE . 0 . )  THK( I , J)=l . 0  •  FOR  NEG.  RHS 


C  UNITS  OF  THICKNESS  ARE  HUNDREDS  OF  M  FOR  CONVENIENCE 
C  SET  INITIAL  WINDS  BEFORE  ALTERATIONS 
U1 ( I , J ) =UN( I , J ) 

VI ( I , J ) =VN ( I , J ) 

C  WEIGHT  WINDS  WITH  THICKNESS  OF  LAYER 
UN( I , J )=U1 ( I , J ) *THK( I , J) 

VN ( I , J ) =V1 ( I , J) *THK ( I , J ) 

40  CONTINUE 

C  IF ( DBUG ( 1 ) )  PRINT  9520 

C  DO  45  J=1 , NROW 

JP=NROW+l-J 

C  45  I F ( DBUG ( 1 ) )  PRINT  9102,  ( THK ( I , JP ) , I  =  1 1 , 1 2 ) 

C  COMPUTE  ORIGINAL  DIVERGENCE  AND  VORTICITY 

DO  170  J=1 , NROWM1 
DO  170  1=1 , NCOLM1 
UE=0 . 5* ( UN ( I + 1 , J ) +UN ( I + 1 , J+l ) ) 

UW=0 . 5* ( UN ( I , J )  +UN ( I , J 1 1 ) ) 

VSO=0. 5»(VN{ I+1,J)+VN( I,J) ) 

VNO=0 . 5* ( VN ( I , J+l ) +VN { I  +  l ,  J  +  l ) ) 

DUE=GSI * { UE-UW ) 

DVN=GSI* ( VNO-VSO) 

DI ( I , J ) =DUE+DVN  •  DIV,  UNITS  ARE  10-6  SEC-1 
IF  { IVORT  .EQ.  2)  GO  TO  170 
VE=0. 5*(VN( 1+1, J)+VN( 1  +  1, J+l) ) 

VW=0 . 5* ( VN( I , J )  + VN ( I , J  + 1 ) ) 


r.on  no 


USO=0 . 5  * ( UN ( I  + 1 , J ) +UN ( I , J ) ) 

UNO=0. 5*(UN< I,J+1)+UN( 1+1, J+l) ) 

DVE=GSI* (VE-VW) 

DUN=GSI* (UNO-USO) 

VT( I , J)=DVE-DUN  !  VORTICITY 

1 70  CONTINUE 

I F ( DBUG ( 1 ) )  PRINT  9570 
DO  173  J= 1 , NROW 
JP=NROW+ 1“J 

173  I F ( DBUG { 1 )  )  PRINT  9100,  ( VT( I , JP ) , 1= 1 1 , 12 ) 

IF ( DBUG ( 1 } )  PRINT  9571 
DO  172  J=1 , NROW 
JP=NROW+l-J 

C  172  IF( DBUG ( 1 ) )  PRINT  9100,  (Dl( I, JP) , 1=11, 12) 

C  FOR  NONDIVERGENCE  SET  DDIJ=0.0 
DDI J=0 . 0 
LG=0 

RA=0 . 4  !  RELAXATION  FACTOR 

210  LG=LG+ 1 

DO  500  J=1 , NROWMl 
DO  500  1=1 , NCOLM1 

UE=0 . 5* ( UN ( I +1 , J ) +UN ( 1+1 , J+l ) ) 

UW=0. 5«(UN(I,J)  +UN( I , J+l ) ) 

VSO=0. 5*(VN( I +1 , J ) +VN ( I , J ) ) 

VNO=0. 5*(VN{ I,J+1)+VN(I+1, J+l) ) 

DUE=GS1*  (UE-UW) 

DVN=GSI * ( VNO-VSO ) 

DI  ( I , J ) =DUE+DVN 

CUIJ=.05*GS  * ( DDI J-DI ( I , J ) ) *RA 
CVIJ= . 05*GS* ( DDIJ-DI { I , J ) ) *RA 

IF  (CUIJ  .LT.-1.0)  CUIJ=-1.0  !  LIMIT  CHANGES 

IF  (CUIJ.GT.  1.0)  CUIJ=1.0 

IF  (CVIJ  . LT .-1.0)  CVIJ=-1.0 

IF  (CVIJ.GT.  1.0)  CVIJ=1.0 

UN ( 1+1 , J ) =UN ( I +1 , J ) +CUIJ 

UN ( I  + 1 , J  + 1 ) =  UN ( I  + 1 , J  + 1 )  +CUIJ 

UN ( I , J ) =UN ( I , J )  -CUIJ 

UN  ( I , J+ 1 ) =UN ( I , J  + 1 )  -CUIJ 

VN  (  1  +  1,  J)=VN(  I.+  l,  J)-CVIJ 

VN ( I , J ) =VN ( I , J ) -CVIJ 

VN ( I , J+l ) =VN ( I , J+l ) +CVIJ 

VN( 1+1 , J+l ) =VN ( 1+1 , J+l ) +CVI J 

IF  (IVORT  .EQ.  2)  GO  TO  490 

VE=0 . 5* ( VN ( 1+1 , J ) +VN ( 1+1 , J+l ) ) 

VW=0 . 5* ( VN ( I , J ) +VN ( I , J+l ) ) 

USO=0. 5*(UN(I+1,J)+UN{I,J) ) 

UNO=0 . 5* ( UN ( I , J+ 1 ) +UN ( 1+ 1 , J+l ) ) 

DVE=GSI* (VE-VW) 

DUN=GSI * ( UNO-USO) 

VO( I , J ) =DVE-DUN 

CVIJ= . 05  +  GS  * ( VT ( I , J ) -VO ( I , J ) ) * RA 
CUIJ=.05‘GS*(VT( I , J ) -VO( I , J ) ) *RA 
IF  (CUIJ  . LT .-1.0)  CUIJ=-1.0  •  LIMIT  CHANGES 

IF  (CUIJ.GT.  1.0)  CUIJ=1.0 
IF  (CVIJ  .LT.-1.0)  CVIJ=-1.0 
IF  (CVIJ.GT.  1.0)  CVIJ=1.0 
UN  ( I , J ) =UN ( I , J )  +CUIJ 
UN( 1+1, J )=UN( t+1 , J ) +CUIJ 
UN( I , J+l ) =UN( I, J+l) -CUIJ 
UN( 1+1, J+1)=UN( 1+1, J+l )-CUIJ 
VN ( 1+1, J)=VN( I+1,J)+CVIJ 
VN ( I + 1 , J+ 1 ) =VN ( I + 1 , J+l )  +CV1J 
VN ( I , J ) =VN ( 1 , J ) -CV I J 
VN ( I , J+ 1 ) =VN ( I , J+l )-CVIJ 
490  CONTINUE 

C  TO  KEEP  WINDS  0.0  WHERE  RHS  IS  NEGATIVE 
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IF  ( RHS ( I , J , L ) . LE .0.0)  THEN 
UN ( I , J ) =0 . 0 
VN ( I , J )=0 . 0 
END  IF 

C  HOLD  ANALYZED  WIND  COMPONENTS  FIXED  AT  PTS  WHERE  IFXPT=1 
IF  (IFXPT(I,J)  .NE.  1)  GO  TO  500 
UN (I,J)=U1(I,J) *THK( I , J ) 

VN  ( I ,  J  ) =V1 ( I , J ) *THK ( I ,  J  ) 

500  CONTINUE 

IF  (LG. GT. NITER)  GO  TO  540 
DO  530  J=2,NROWMl 
DO  530  I=2,NCOLMl 

IF  (ABS(DDIJ-DI ( I , J ) ) .GT .50.0)  GO  TO  210 
C  IF  (ABS{ VT( I , J ) -VO( I ,  J ) ) .GT .50.0)  GO  TO  210 

530  CONTINUE 
540  CONTINUE 

C  I F ( DBUG ( 1 ) )  PRINT  9570 

C  IF( DBUG ( 1 } )  PRINT  9580 

C  DO  510  J=1 , NROW 

JP=NROW+l-J 

C  510  IF( DBUG ( 1 ) )  PRINT  9100,  ( VO( I , JP ) , 1  =  11 , 1 2 ) 

C  IF( DBUG ( 1 ) )  PRINT  9200,  LG 

C  IF( DBUG( 1 ) )  PRINT  9571 

C  IF ( DBUG ( 1 ) )  PRINT  9580 

C  DO  520  J=1 , NROW 

jp-  NROW't’l-  J 

C  520  IF(DBUG( 1 ) )  PRINT  9100,  ( DI ( I , JP ) , 1  =  II , 1 2 ) 

SUM1=0.0 
SUM2=0 . 0 
01=0.0 

DO  1040  J=1 , NROW 
DO  1040  1=1 , NCOL 

IF  ( RHS ( I , J , L ) . LE .0.0)  GO  TO  1040  !  OMIT  THESE  PTS 

01=01+1.0 

UN ( I , J ) =UN ( I , J ) /THK ( I , J ) 

VN ( I , J ) =VN ( I ,  J ) /THK ( I , J ) 

U1(I,J)=U1(I,J) -UN( I , J ) 

VI  ( I ,  J )  =V1  ( I ,  J )  -VN  ( I ,  J ) 

SUM1  =  SUM1+U1(.I,J) 

SUM2=SUM2+V1 ( I , J ) 

1040  CONTINUE 

SUM1 = SUM1/Q 1 
SUM2=SUM2/Q1 

C  NORMALIZE  ORIG.  AVERAGE  VALUES 
C  DO  1045  J=1 , NROW 

C  DO  1045  1=1, NCOL 

C  UN ( I , J ) =UN ( I , J ) +SUM1 

C  VN ( I , J ) =VN ( I , J )  +SUM2 

IF(DBUG( 1 ) )  PRINT  1145,  L 
DO  1160  J=l,NROW 
JP=NROW+l-J 

1160  IF(DBUG( 1 ) )  PRINT  9150,  ( Ul( I , JP ) , 1=11 , 12 ) 

IF(DBUG( 1 ) )  PRINT  1155 
DO  1170  J=1 , NROW 
JP=NROW+l-J 

1170  I F ( DBUG ( 1 ) )  PRINT  9150,  ( VI ( I , JP) , 1  =  11 , 1 2 ) 

C  WRITE  LEVEL  4  DIV  WINDS  TO  OUTPUT  FILE 
IF  (L  .NE.  4)  GO  TO  700 
IF( . NOT . RYT3 )  GO  TO  700 
DO  690  J=l,  NROW 
JR=NROW+l-J 

690  WRITE  (3,9065)  ( U1  (  I , JR ) , I  =  1 , NCOL) 

DO  695  J=l,  NROW 
JR=NROW+l-J 

695  WRITE  (3,9065)  ( VI ( I , JR ) , 1=1 , NCOL) 

700  CONTINUE 
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C  CHANGE  BACK  TO  3D  ARRAYS 

C 

DO  580  J  =  1 , NROW 
DO  580  1  =  1 , NCOL 

U( I , J, L) =UN( I / J ) 

V  ( I ,  J ,  L ) =VN ( I ,  J ) 

580  CONTINUE 
600  CONTINUE 

IF(DBUG( 1 ) )  PRINT  9003 

1145  FORMAT ( ’  U  COMP.  DIVERGENT,  LEVEL  =’I3) 

1155  FORMAT ( ’  V  COMP.  DIVERGENT  ') 

9002  FORMAT  (//'  BEGIN  SUBROUTINE  BAL5  ’/) 
yu03  FORMAT  (//’  END  OF  SUBROUTINE  BAL5  '/) 

9065  FORMAT  (10F8.1) 

9150  FORMAT ( '  ’ ,33F4.1) 

9562  FORMAT  (/'  POINT  WITH  FIXED  WIND  IS  COL  ='I3,'ROW  ='I3/) 

RETURN 

END 

C 

C**t***«*************************************************************** 

c 

SUBROUTINE  DOPSIG ( DBUG , KEY ) 

C  ASSIGN  DOPPLER  WIND  PROFILES  TO  SIGMA  SURFACES. 

C  MISSING  WINDS  ARE  DENOTED  BY  -999. 

C  IF  SOUNDING  IS  NOT  COMPLETE  THE  LAST  REPORTED  WIND 
C  IS  USED  AT  THE  HIGHEST  ALTITUDES. 

C  WIND  DATA  START  AT  40  M  AND  CONTINUE  AT  30-M  INTERVALS 
C  TO  610  M.  THERE  ARE  20  POINTS  IN  A  PROFILE. 

C  FILL  IN  MISSING  DATA  WITH  NEAREST  POINT  UP  OR  DOWN. 

C  FOR  PROGRAM  DIABWND,  JAN  '85. 

C  BY  R.M.  ENDLICH,  SRI  INTN’L,  MENLO  PARK  CA  94025. 

C  MODIFIED  BY  FLLUDWIG— 12/85 

LOGICAL  KEY , DBUG( 15 ) , UTMUSE 

DIMENSION  DPHT( 5,50),  DPUC(5,50),  DPVC( 5 , 50 ) , RHS1 ( 5 , 6 ) 

DIMENSION  DPWD(5,50) , DPWS( 5 , 50 ) , XS ( 5 ) ,YS(5) ,STLT(5) ,STLN(5) 

COMMON  /STALOC/  XG(50),  YG(50) 

COMMON  /WINDS/  USIG(50,6),  VSIG(50,6) 

COMMON  /LIMITS/NCOL, NROW, NLVL, NCOLMl , NROWM1 , 

$  LOWIX( 5,3), LOWIY (5,3), SFCMAX 

COMMON/RARS/RHS (25,25,6) 

COMMON/PARMS/ZTOP , DS , DSIGMA , NLVLM1 , XHT1 , XHT2 , XI , Y1 , 

1  X2,Y2,UG,VG, RATIO, TDSI 
COMMON  /SITE/  IXS,  JYS,  THSITE,  IGRID 
COMMON  /ANCHOR/  SLAT,  SLNG  , UTMUSE 
COMMON  /NUMOBS/  NUMDOP ,  NUMNWS,  NUMTOT 
COMMON  /UPWIND/  UTOP ,  VTOP 
INTEGER  STAID( 5 ) 

C  VARIABLES  ARE: 

DPUC=U  COMPONENT  OF  LajPPLEK  WIND  iN  1U  „ 

C  DPVC=V  COMPONENT  OF  DOPPLER  WIND  IN  MPS 
C  NHTS=NUMBER  OF  POINTS  IN  VERTICAL  WIND  PROFILE 
C  NLVL=NUMBER  OF  SIGMA  LEVELS 
C  RHS=HT  OF  SIGMA  SURFACES  ABOVE  TERRAIN  (M) 

XG , YG=STA.  DIST  IN  X,Y  IN  GRID  UNITS  FROM  SW  CORNER 
IF( DBUG{ 2 ) )  PRINT  9001 
C  FILL  IN  HEIGHTS  OF  DOPPLER  DATA  POINTS 
C - 

C  PRINT  LAT , LONG  OF  ANCHOR  POINT  (UTMUSE=F)  OR  UTM  COORDS  (UTMUSE=T) 

I F ( DBUG ( 2 ) )  PRINT  9007,  SLAT,  SLNG 

C»  READ  LAT  ,  LONG  OF  STATIONS  (DEG)  AND  CONVERT  TO  XS,YS  ( KM) --UTMUSE=F-- 
C  OR  READ  IN  UTM  COORDINATES  (KM)  DIRECTLY--UTMUSE=TRUE--FROM  SW  CORNER 

IF  (KEY)  THEN 

C  READ(12,9014)  NUMDOP 
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IT=1 , NUMDOP ) 


READ(12,»)  NUMDOP 

IF(DBUG( 2) )  PRINT  9027,  NUMDOP 

READ( 12 , * )  (STAID( IT) ,STLT( IT) ,  STLN( IT) , IT=1 , NUMDOP ) 

I F ( DBUG ( 2 ) )  PRINT  9004 
I F ( DBUG ( 2 ) )  PRINT  9014, 

$  (STAID( IT) ,STLT( IT) ,  STLN ( IT ), IT=1 , NUMDOP ) 

DO  15  IT=1, NUMDOP 
IF  (UTMUSE)  THEN 

XS( IT) =STLN { IT ) -SLNG 
YS( IT ) =STLT( IT) -SLAT 
ELSE 

XS( IT) = ( STLN{ IT) -SLNG ) * ( 111 . 2*COS( SLAT/57 . 295) ) 

YS( IT ) = ( STLT( IT) -SLAT) *111 . 2 
END  IF 

XS( IT)=XS( IT) +IXS* ( DS* . 001 ) 

YS( IT) =YS( IT) +JYS* ( DS* . 001 ) 

XG ( IT) =XS( IT)/(DS*.001) 

YG( IT) =YS( IT)/(DS« .001 ) 

15  CONTINUE 

I F ( DBUG ( 2 ) )  PRINT  9008 
I F ( DBUG ( 2 ) )  PRINT  9011, 

$  (XS( IT) ,  YS{ IT) , XG( IT) , YG( IT) , IT=1 , NUMDOP ) 

END  IF 

1  ASSIGN  HTS  OF  SIGMA  SFC  FOR  EACH  SOUNDING 
DO  25  IT=1, NUMDOP 
DO  20  K=1,NLVL 
IX=JNINT(XG( IT) ) 

JY=JNINT( YG( IT) ) 

IF  (IX.LT.l)  IX=1 
IF  (IX.GT.NCOL)  IX=NCOL 
IF  (JY.LT.l)  JY=1 
IF  (JY.GT.NROW)  JY=NROW 
RHS1 ( IT, K ) =RHS( IX, JY , K ) 

20  CONTINUE 

I F ( DBUG ( 2 ) )  PRINT  9024 

IF(DBUG(2))  PRINT  9015,  IT,  (RHS1( IT,K) ,K=1,NLVL) 

25  CONTINUE 

C  READ  DOPPLER  SOUNDINGS  FOR  NUMDOP  STATIONS 
DO  450  IT=1, NUMDOP 
IF  (KEY)  THEN 

C  READ  (12,9014)  NHTS 

READ  (12,*)  NHTS 

I F ( DBUG ( 2 ) )  PRINT  9003,  IT,  NHTS 

C  READ  (12,9009)  ( DPHT( IT, LL) , DPWD( IT, LL ) ,  DPWS(IT,LL) 

C  $  LL=1 , NHTS ) 

READ  (12,*)  ( DPHT{  IT ,  LL )  ,  DPWD(  IT ,  LL )  ,  DPWS(IT,LI,), 

$  LL- 1 , NUTS j 

I F ( DBUG ( 2 ) )  PRINT  9009, 

$  ( DPHT( IT, LL) , DPWD( IT, LL ) ,  DPWS( IT, LL) , LL=1 , NHTS ) 

C  CHANGE  DIRECTION  AND  SPEED  (MPS)  TO  U  AND  V 
DO  40  LL=1 , NHTS 

IF  ( DPWD( IT, LL)  .EQ. -999.0)  THEN 
DPUC( IT, LL)=-999 . 0 
DPVC ( IT, LL ) =-9  9  9 . 0 

ELSE 

DPUC( IT, LL) =-DPWS( IT, LL ) *SIN ( DPWD( IT , LL ) /57 . 295 ) 
DPVC ( IT, LL) =-DPWS ( IT, LL ) *COS ( DPWD( IT , LL ) /57 . 29 5 ) 

END  IF 

40  CONTINUE 

I F ( DBUG ( 2 ) )  PRINT  901 0 , ( DPUC( IT, LL) , 

$  DPVC ( IT, LL ) , LL= 1 , NHTS ) 

C  FILL  IN  MISSING  DATA  IF  NEEDED 
C  START  FROM  BOTTOM  AND  GO  UPWARD 
LL=0 

60  LL=LL  +1 

IF  (LL  .EQ.  NHTS)  GO  TO  100 


( RHS1 ( IT, K ) , K=1 , NLVL ) 
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IF  { DPUC ( IT, LL)  .EQ. -999.0)  THEN 
GO  TO  70 
ELSE 

GO  TO  60 
END  IF 
70  LL1=LL 

75  LL1=LL1+1 

IF  l DPUC ( IT, LL1 )  .EQ. -999.0  .AND.  LL1  .EQ.  NHTS ) 

+  GO  TO  100 

IF  (DPUC( IT, LL1)  .EQ. -999.0)  THEN 
GO  TO  75 
ELSE 

DPUC ( IT, LL ) =DPUC ( IT,LL1) 

DPVC ( IT, LL)=DPVC( IT, LL1 ) 

IF( DBUG( 2 ) )  PRINT  9015,  LL,  DPUC(IT,LL),  DPVC(IT,LL) 

GO  TO  60 
END  IF 

100  CONTINUE 

I F ( DBUG ( 2 ) )  PRINT  9010 , ( DPUC( IT, LL ) , DPVC( IT, LL ) , LL=1 , NHTS ) 
C  START  FROM  TOP  AND  GO  DOWN 
LL=NHTS+1 
110  LL=LL-1 

IF  (LL  -EQ.  1)  GO  TO  140 
IF  (DPUC( IT, LL)  .EQ. -999.0)  THEN 
GO  TO  120 
ELSE 

GO  TO  110 
END  IF 
120  LL1=LL 

125  LL1=LL1-1 

IF  (DPUC( IT, LL1 )  .EQ. -999.0)  THEN 
GO  TO  125 
ELSE 

DPUC ( IT, LL) =DPUC ( IT , LL1 ) 

DPVC ( IT, LL) =DPVC ( IT, LL1 ) 

I F ( DBUG ( 2 ) )  PRINT  9015,  LL,  DPUC(IT,LL),  DPVC(IT,LL) 

GO  TO  110 
END  IF 
140  CONTINUE 

IF( DBUG ( 2 ) )  PRINT  9012 

I F ( DBUG ( 2 ) )  PRINT  9010 ,( DPUC< IT, LL ), DPVC ( IT, LL ), LL=1 , NHTS ) 
END  IF 

BEGIN  INTERPOLATION  SCHEME 

DO  400  K=2 ,  NLVL  !  COUNTER  FOR  SIGMA  LEVELS 
LL=0 

200  LL=LL+1 

IF  ( RHS1 ( IT, K ) . GT .0.0)  GO  TO  365 
USIG ( IT, K ) =0 . 0 
VSIG( IT, K ) =0 . 0 
GO  TO  400 

365  IF  (RHS1 ( IT, K ) . LE . DPHT( IT, LL) )  GO  TO  320 

IF  (RHS1( IT,K) ,GE.DPHT( IT, NHTS) )  GO  TO  380 
IF  (RHS1 ( IT, K ) .GE . DPHT( IT, LL)  .AND.  RHS1 ( IT, K ) . LE . 

+  DPHT( IT, LL+1 ) )  GO  TO  360 

GO  TO  200 

('  FOR  LEVELS  BELOW  1ST  MEASURED  WINDS  (ASSUME  SPEED=0  AT  1M 

C  AND  AT  GROUND  WHERE  RHS=0) 

320  USIG( IT, K)=(DPUC( IT, LL) ) * ( ALOG10 ( RHS1 ( IT, K ) ) )/ 

+  ALOG10(DPHT( IT,  LL) ) 

VSIG( IT, K)  =  (DPVC( IT, LL) ) * ( ALOG10 ( RHSl ( IT,  K ) ) )/ 

+  ALOG10(DPHT( IT, LL) ) 

GO  TO  400 

C  FOR  SIGMA  LEVELS  BETWEEN  DOPPLER  WIND  POINTS 

360  RATIO= ( ALOG10 ( RHSl ( IT , K ) ) -ALOG1 0 ( DPHT( IT, LL) ) )/ 

+  ( ALOG10 ( DPHT( IT, LL+ 1 ) ) -  ALOGIO ( DPHT( IT, LL) ) ) 

USIG( IT , K ) =DPUC ( IT , LL )  + ( DPUC( IT , LL+ 1 ) -DPUC ( IT , LL) ) 
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+  ‘RATIO 

VSIG( IT, K)=DPVC( IT, LL)  + (DPVC( IT, LL+1 ) -DPVC ( IT, LL) ) 

+  ‘RATIO 

GO  TO  400 

C  FOR  SIGMA  LEVELS  ABOVE  LAST  DOPPLER  POINT 
380  USIG( IT, K)=DPUC( IT,NHTS) 

VSIG ( IT, K ) =DPVC ( IT, NHTS ) 

400  CONTINUE 

IF( DBUG( 2 )  )  PRINT  9020,  IT 

IF(DBUG{ 2 ) )  PRINT  9010,  (USIG(IT,K),  VSIG(IT,K),  K-1,NLVL) 
450  CONTINUE 

C  GET  AVER  UPPER  WIND  FOR  POSSIBLE  USE  IN  GEOSIG  IF  GEOS  WIND 
C  IS  NOT  AVAILABLE 
UTOP=0 
VTOP=0 
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9001 

9u02 

9003 

9004 

9007 

9008 

9009 

9010 

9011 

9012 

9014 

9015 
9020 
9024 

9026 

9027 


DO  460  IT-1,  NUMDOP 

UTOP=UTOP+USIG ( IT , NLVL ) 

VTOP=VTOP+VSIG( IT,NLVL) 

CONTINUE 

UTOP=UTOP/ ( FLOAT ( NUMDOP ) ) 

VTOP=VTOP/( FLOAT (NUMDOP ) ) 

I F { DBUG ( 2 ) )  PRINT  9026,  UTOP,  VTOP 
I F ( DBUG ( 2 ) )  PRINT  9002 

BEGIN  SUBROUTINE  DOPSIG  ’/) 

END  OF  SUBROUTINE  DOPSIG  ’/) 

STA.  =  1 13 , 1  NO.  OF  POINTS  IN  PROFILE= ’ 13/ ) 
LATITUDE  AND  LONGITUDE  OF  STATIONS’/) 

THE  ANCHOR  PT .  IS  AT  LAT  = ’ F9 . 3 , ’  AND  LONG= ’ F9 . 3 ) 
X  AND  Y  OF  STATIONS  IN  KM  AND  IN  COLS, ROWS  FROM  SW 
CORNER  OF  COARSE  GRID’/) 

FORMAT  ( 3X, 3F7 . 1 ) 

FORMAT  ( 3X, 8F6 . 1 ) 

FORMAT  (/4X,2F11.3, 8X, 2F10 . 2 ) 

FORMAT  (/’  SOUNDING  WITH  FILLED  IN  DATA  ’/) 

FORMAT  ( 3X, 15 , 2F8 . 2 ) 

FORMAT  ( 3X, 15 , 6F8 . 2 ) 

FORMAT  (/’  STATION  NUMBER  =’I6/) 

FORMAT  {/’  HEIGHTS  OF  SIGMA  SURFACES  ’/) 

FORMAT  (/’  AVER.  UPPER  WIND.  U  =  ’  F6  1.’  V= ’ F6 . 1 / ) 

FORMAT  (/’  NUMBER  OF  WIND  SOUNDINGS  - ’ U/ ) 

RETURN 

END 


FORMAT  (/’ 
FORMAT  (/’ 
FORMAT  (/’ 
FORMAT  (//’ 
FORMAT  (//’ 
FORMAT  (//’ 


SUBROUTINE  GEOSIG (DBUG ) 

PREPARE  NWS  HOURLY  REPORTS  OF  WIND  DIRECTION,  WIND  SPEED 
(KNOTS),  SEA  LEVEL  PRESSURE  (MB),  AND  TEMPERATURE  (DEG  F) 
FOR  INPUT  TO  WIND  ANALYSIS  FOR  DIABLO  PGE  SITE. 

COMPUTE  GEOS  WIND  FROM  PRESSURE  AT  THREE  STATIONS  AND 
CORRECT  IT  FOR  THERMAL  WIND  COMPONENT  (IF  DESIRED). 

ASSUME  THAT  WIND  COMPONENTS  VARY  WITH  LOG (HEIGHT)  BETWEEN 
ANEMOMETER  HT  AND  GEOS  WIND  AT  HT  GEOSHT( ABOUT  500M) . 

AT  EACH  NWS  STATION  INTERPOLATE  WINDS  TO  SIGMA  SURFACES. 
BY  RM  ENDLICH,  SRI  INTN’L,  MENLO  PARK  CA  94025  DEC  ’84. 
VAR IABLES 

NUMNWS  =  NUMBER  OF  NWS  REPORTS 
NWSID  =  IDENTIFICATION  ID  OF  NWS  STATION 
WD  =  WIND  DIRECTION 
SP  =  WIND  SPEED 

STLT  =  STATION  LATITUDE  IN  DEGS  AND  HUNDREDTHS 
STLN  =  STATION  LONGITUDE 

s.  PRESS  =  STATION  SEA  LEVEL  PRESSURE  IN  INCHES  HG 
C  TEMP  =  STATION  TEMP  IN  DEG  F 
LOGICAL  DBUG (15),  UTMUSE 
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INTEGER  STAID( 50 ) 

DIMENSION  STLT( 50 ) ,  STLN(50),  PRESS (50) 

DIMENSION  TEMP( 50 ) 

DIMENSION  RHS1( 50, 6),  UNWS(50,6),  VNWS(50,6) 

DIMENSION  UCOMP (50),  VCOMP(50),  WD(50),  WS(50) 

DIMENSION  XS( 50 ) ,  YS(50) 

DIMENSION  XG1{ 50),  YG1(50),  USIG1(50,6),  VSIG1(50,6) 

COMMON  /LIMITS/NCOL, NROW, NLVL , NCOLM1 , NROWM1 , 

$  DOWIX( 5,3), LOWI Y (5,3), SFCMAX 

COMMON/RARS/RHS ( 25,25,6) 

COMMON/PARMS/ZTOP , DS , DSIGMA, NLVLM1 , XHT1 , XHT2 , XI , Y1 , 

1  X2,Y2,UG,VG, RATIO, TDSI 

COMMON  /SITE/  IXS,  JYS,  THSITE,  IGRID 
COMMON  /STALOC/  XG ( 50 ) ,  YG ( 50 ) 

COMMON  /WINDS/  USIG(50,6),  VSIG(50,6) 

COMMON  /ANCHOR/  SLAT,  SLNG  ,UTMUSE 
COMMON  /NUMOBS/  NUMDOP,  NUMNWS,  NUMTOT 
COMMON  /UPWIND/  UTOP,  VTOP 
I F ( DBUG ( 3 ) )  PRINT  9001 

READ  INPUT  OF  HOURLY  DATA:  STATION  ID,  LATITUDE,  LONGITUDE  (UTMUSE=F)  OR 
UTM  COORDS  ( KM, UTMUSE=T)  PRESSURE,  TEMPERATURE,  WIND  DIR.,  WIND  SPEED  (MPS) 

DO  70  IT=1, NUMNWS 

READ  (12,  10)  STAID(IT),  STLT(IT),  STLN(IT), 

+  PRESS ( IT ) ,  TEMP(IT),  WD( IT) ,  WS(IT) 

READ  (12,*)  STAID ( IT) ,  STLT(IT),  STLN(IT), 

+  PRESS ( IT ) ,  TEMP ( IT ) ,  WD(IT),  WS(IT) 

I F ( DBUG ( 3 ) )  PRINT  20,  IT,  STAID( IT) , STLT( IT) , STLN( IT) , 

+  PRESS ( IT) , TEMP ( IT) 

IF( DBUG ( 3 ) )  PRINT  25,  WD(IT),  WS( IT) 

70  CONTINUE 

CONVERT  LAT , LONG  (UTMUSE=F)  OR  UTM  COORDS ( KM , UTMUSE=T )  TO  XS,YS  (KM) 

0  MEASURED  FROM  SW  CORNER  OF  COARSE  GRID 

DO  15  J=l, NUMNWS 
IF  (UTMUSE)  THEN 

XS ( J ) =STLN ( J ) -SLNG 
YS( J ) =STLT( J ) -SLAT 
ELSE 

XS ( J ) = ( STLN ( J ) -SLNG ) * ( III . 0*COS( SLAT/57 . 295 ) ) 

YS ( J ) = ( STLT( J ) -SLAT) *1 i 1 . 0 
END  IF 

XS(J)=XS(J)  +  IXS* ( DS* . 001 ) 

YS(J)=YS(J)+JYS*(DS* .001) 

XG1 ( J ) =XS ( J )/( DS* . 001) 

YG1(J)=YS(J ) /( DS* . 001) 

15  CONTINUE 

I F ( DBUG ( 3 ) )  PRINT  9008 

IF(DBUG( 3 ) )  PRINT  9011,  (XS(J),  YS ( J ) , XG1 ( J ) , YG1 ( J ) , 

$  J=l, NUMNWS) 

C  CHANGE  DIRECTION  AND  SPEED  (MPS)  TO  U  AND  V 
DO  80  IT=1, NUMNWS 

UCOMP ( IT ) =-WS( IT) *SIN(WD( IT)/57 . 295) 

VCOMP ( IT)=-WS( IT) *COS(WD( ITJ/57.295) 

I F ( DBUG ( 3 ) )  PRINT  9002 

I F ( DBUG ( 3 ) )  PRINT  9013,  IT,  UCOMP(IT),  VCOMP(IT) 

80  CONTINUE 

C  ASSIGN  HEIGHTS  TO  SIGMA  SURFACES 
DO  125  J=l, NUMNWS 
DO  120  K= 1 , NLVL 
IX=JNINT  (XG1(J) ) 

JY=JNINT  (YG1(J)) 

IF  (IX.LT.l)  IX=1 
IF  (IX.GT.NCOL)  IX=NCOL 


60 


IP  (JY.LT.l)  JY=1 
IF  (JY.GT.NROW)  JY=NROW 
RHS1 ( J ,  K ) =RHS{ IX , JY , K ) 

120  CONTINUE 

IF(DBUG( 3 ) )  PRINT  9024 

I F ( DBUG ( 3 ) )  PRINT  9015,  J,  ( RHS1 ( J , K ) , K=1 , NLVL) 

125  CONTINUE 

C  IF  LESS  THAN  3  NWS  STATIONS  CANT  COMPUTE  GEOS  WIND.  INSTEAD 
C  USE  UTOP,  VTOP  FROM  DOPPLER. 

IF  ( NUMNWS . LT . 3 )  UGEOS=UTOP 
IF  ( NUMNWS . LT . 3 )  VGEOS=VTOP 
IF  (NUMNWS.LT. 3)  GO  TO  200 
C 

C* ‘THIS  SECTION  COMPUTES  GEOSTROPHIC  WINDS  AND  IT  ALSO  CORRECTS 
C  THEM  FOR  THERMAL  WINDS  IF  NTHERM=1 . 

C 

ISl=NUMNWS-2 
IS2=NUMNWS-1 
IS3=NUMNWS 
C  SET  NTHERM 

NTHERM= 1 

PR1=PRESS( ISl ) * ( 1013 . 3/29 . 92 ) 

PR2=PRESS( IS2 ) * ( 1013 . 3/29 . 92 ) 

PR3=PRESS( IS3 ) *( 1013 . 3/29 . 92 ) 

STLT1=STLT( ISl) 

STLT2=STLT( IS2 ) 

STLT3=STLT( IS3 ) 

STLN1=STLN ( ISl ) 

STLN2=STLN ( IS2 ) 

STLN3=STLN ( IS3 ) 

TMPl=TEMP ( ISl ) 

TMP2=TEMP(IS2) 

TMP3=TEMP( IS3 ) 

C  DEFINE  CONSTANTS 

AV LAT= .333*  ( STLT 1 + STLT  2 + STLT3 ) / 5  7 . 2  9  5 
FC=14 . 584 *SIN  (AVLAT) 

C  CORIOLIS  FORCE  IN  UNITS  10-5  SEC-1 
COSLAT=COS  (AVLAT) 

DENOM= ( STLT2-STJ/T1 ) * ( STLN3-STLN1 ) - ( STLT3-STLT1 )  * 

+  (STLN2-STLN1) 

RHO=l . 1 

C  DENSITY  IN  UNITS  10-3  G/CM3 
C2=100 . 0/(RHO*1 . 112 ) 

C  COMPUTE  GEOSTROPHIC  WINDS 

DPDLT= ( ( STLN2-STLN1 ) * ( PR3-PR1 )- ( STLN3-STLN1 )  * 
f  ( PR2-PR1 ) )/( -DENOM) 

DPDLN= ( (STLT2-STLT1 ) * ( PR3-PR1 ) - ( STLT3-STLT1 )  * 

+  ( PR2-PR1 ) ) /DENOM 

UGEOS= ( -C2/FC ) *DPDLT 
VGEOS= (C2/FC ) * ( DPDLN/COSLAT ) 

C  SPEED  UNITS  ARE  M  SEC-1 

I F ( DBUG ( 3 ) )  PRINT  30,  UGEOS,  VGEOS 
C  THIS  PART  MAKES  THERMAL  WIND  CORRECTION  TO  UGEOS,  VGEOS 
C  NTHERM  IS  INDICATOR  FOR  USE  (WHEN=1) 

IF  (NTHERM  .NE.  1)  GO  TO  200 

FVNIN=5 . 0/9 . 0 

TMP1= (TMP1-32 . 0 ) *FVNIN 

TMP2= ( TMP2-3  2 . 0 ) *FVNIN 

TMP3= ( TMP3 - 3 2 . 0 ) *FVNIN 

AVTMP=273 . 0+ . 333* (TMP1+TMP2+TMP3 ) 

0  GRAVITY=9.8  M  SEC- 2,  TEMP  IN  DEG  K 
C3=9.8/(FC*1.112) 

DTDLT= ( ( STLN2-STLN1 ) *  (TMP3-TMP1 ) - ( STLN3-STLN1 ) 

4  *(TMP2-TMP1) )/(-DENOM) 

DTDLN= ( ( STLT2-STLT1 ) *  (TMP3-TMP1 ) - ( STLT3-STLT1 ) 

+  *(TMP2-TMP1) ) /DENOM 
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UTHERM=- ( C3/AVTMP ) *DTDLT 
VTHERM= ( C3/AVTMP ) » ( DTDLN/COSLAT ) 

C  ASSUME  SHEAR  ACTS  OVER  LAYER  OF  DEPTH  DEPL  (IN  M) 

DEPL=200  !  TIE  THIS  IN  TO  AVTHK  FOR  B  LYR  TOP 
USHEAR=UTHERM*DEPL 
VSHEAR  = VTHERM  *  DEPL 

I F ( DBUG ( 3 ) )  PRINT  40,  USHEAR ,  VSHEAR,  DEPL 
UGEOS=UGEOS+USHEAR 
VGEOS=VGEOS+ VSHEAR 
I F ( DBUG ( 3 ) )  PRINT  50 
IF( DBUG ( 3 ) )  PRINT  30,  UGEOS,  VGEOS 
200  CONTINUE 

*  *END  OF  GEOSTROPHIC  WIND  SECTION 
INTERPOLATE  WINDS  BETWEEN  SFC  AND  GEOSHT 
DO  450  IT=1 ,  NUMNWS 

DO  400  K=2 ,  NLVL  !  COUNTER  FOR  SIGMA  LEVELS 
IF  (RHSl ( IT, K ) . GT .0.0)  GO  TO  370 
UNWS( IT , K ) =0 . 0 
VNWS ( IT, K ) =0 . 0 
GO  TO  400 

370  IF  ( RHSl ( IT, NLVL) . GE . 800 . )  GEOSHT=800. 

IF  (RHS1( IT, NLVL) .LT. 800. )  GEOSHT=RHSl ( IT, NLVL) 

IF  (RHS1( IT, K) .GT. GEOSHT)  GO  TO  380 
C  FOR  SIGMA  LEVELS  BETWEEN  SURFACE  OBS  AND  GEOSHT 

RATIO= ( ALOG1 0 ( RHSl ( IT , K ))- 1 . 0 )/( ALOG 1 0 ( GEOSHT )- 1 . 0 ) 

UNWS ( IT, K ) =UCOMP( IT)  + ( UGEOS-UCOMP ( IT ) ) ‘RATIO 
VNWS ( IT , K ) =VCOMP ( IT )  + ( VGEOS-VCOMP ( IT ) ) ‘RATIO 
GO  TO  400 

C  FOR  SIGMA  LEVELS  ABOVE  GEOSHT 
360  UNWS ( IT , K ) =UGEOS 

VNWS ( IT, K ) =VGEOS 
400  CONTINUE 

I F ( DBUG ( 3 ) )  PRINT  9020,  IT 

IF( DBUG ( 3 ) )  PRINT  9010,  (UNWS(IT,K),  VNWS(IT,K),  K=1,NLVL) 
450  CONTINUE 

C  INCLUDE  NWS  DATA  WITH  DOPPLER  DATA  IN  ARRAYS  NEEDED 
C  FOR  OBJECTIVE  ANALYSIS. 

DO  460  N=l,  NUMNWS 
M=N+NUMDOP 
XG ( M ) =XG1 ( N ) ' 

YG(M)=YG1(N) 

DO  460  K=l,  NLVL 

USIG (M, K ) =UNWS ( N , K ) 

VSIG(M, K ) =VNWS ( N , K ) 

460  CONTINUE 

NUMTOT=NUMDOP +NUMNWS 
DO  480  M=l,  NUMTOT 
I F ( DBUG ( 3 ) )  PRINT  9020, M 

I F ( DBUG ( 3 ) )  PRINT  9010,  (US1G(M,K),  VSIG(M,K),  K=1,NLVL) 
480  CONTINUE 

IF(DBUG( 3 ) )  PRINT  9003 
10  FORMAT  ( 15 , 2F8 . 2 , 4F7 . 1 ) 

20  FORMAT  (/'  NO.='I3,’  STA  ID='I5,'  LAT=’F7.2,'  LONG. 

+  = ' F7 . 2 , '  PRESS='F7 . 2, '  TEMP='F6.1/) 

25  FORMAT  ('  WIND  DIR='F6.1,'  WIND  SPEED= ' F5 . 1/ ) 

30  FORMAT  ('  GEOSTROPHIC  WIND  COMPONENTS,  MPS,  U  = ’ F5 . 1 , 

+  *  V  = ' F5 . 1/ ) 

40  FORMAT  (//'  THERMAL  WIND  SHEAR  COMPONENTS,  MPS,  U  =' 

+  F5.1,’  V  =’F5.1,'  LAYER  DEPTH  = ' F6 . 1 , '  M'/) 

50  FORMAT  (//'  GEOSTROPHIC  WIND  CORRECTED  FOR  THERMAL  WIND’) 
5001  FORMAT  (/'  BEGIN  SUBROUTINE  GEOSIG  '/) 

0002  FORMAT  (//'  STA.  NO.  U  COMP  V  COMP  '/) 
tOO 3  FORMAT  (/'  END  OF  SUBROUTINE  GEOSIG  '/) 

t008  FORMAT  (//'  X  AND  Y  OF  STATIONS  IN  KM  AND  IN  COLS , ROWS 
t  FROM  SW  CORNER  OF  COARSE  GRID'/) 

-1010  FORMAT  (  3X,  8F6 . 1 ) 
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FORMAT  (/4X,2F11.0,8X,2F10.1) 

FORMAT  (/'  STA.  ='I5,'  UCOMP= ' F8 . 1 , ’  VC 
FORMAT  (3X,I5,8F8.1) 

FORMAT  (/'  STATION  NUMBER  =  *  16/) 

FORMAT  (/'  HEIGHTS  OF  SIGMA  SURFACES  '/) 
RETURN 
END 


VCOMP= ' F8 . 1/ ) 


SUBROUTINE  GPAN(DBUG) 

:  THIS  ROUTINE  MAKES  A  GRID  PT  ANALYSIS  FROM  AVAILABLE 
:  WIND  OBSERVATIONS.  THE  OBSERVATIONS  HAVE  BEEN  INTERPO- 
:  LATED  TO  SIGMA  SURFACES.  THE  WEIGHTING  GIVEN  TO  EACH 
;  STATION  IS  INVERSELY  PROPORTIONAL  TO  ITS  DISTANCE  FROM 
:  THE  GRID  POINT. 

:  TOR  SUBROUTINE  DIABWND,  JAN  '85. 

BY  R.M.  ENDLICH ,  SRI  INTN'L,  MENLO  PARK  CA  94025. 

LOGICAL  DBUG( 15 ) 

INTEGER  STAID{ 6 ) 

DIMENSION  WT( 50 ) 

COMMON  /WINDS/  USIG(50,6),  VSIG(So.f.) 

COMMON  /STALOC/  XG(50),  YG(5u) 

COMMON  /LIMITS/NCOL, NROW , NLVL , NCOLM1 , NROWM1 , 

S  LOWIX( 5,3), LOWI Y (5,3), SFCMAX 

COMMON /UARS/U (25,25,6) , UA( 25,25,6) , V ( 25 , 25 , 6 ) , VA( 25,25,6) 
COMMON  /NUMOBS/  NUMDOP,  NUMNWS ,  NUMTOT 
COMMON  /PRLIM/  II,  12 
C  VARIABLES  ARE: 

C  USIG=U  COMP  OF  WIND  ON  A  SIGMA  SFC 

C  VSIG=V  COMP  OF  WIND  ON  A  SIGMA  SFC 
C  WT= WEIGHT  ASSIGNED  TO  A  GIVEN  STATION 

C  U( I , J , K )  AND  V( I , J , K )  ARE  FINAL  WIND  COMPONENTS 
C  NLVL=NUMBER  OF  SIGMA  LEVELS 

C  XG, YG  ARE  DISTANCES  OF  STATIONS  FROM  SW  CORNER  OF  COARSE 
C  GRID  AND  ARE  MEASURED  IN  GRID  UNITS  (DSCRS) 

C  IN  ARRAYS  (IT,K)  IT  DENOTES  STATIONS,  K  DENOTES  LEVELS 

C  IN  3-D  ARRAYS  I ,  J1,  K  DENOTE  COLS,  ROWS,  LEVELS  FROM  SW  CORNER 


I F ( DBUG ( 4 ) )  PRINT  9001 

SET  MINIMUM  DISTANCE  (GRID  UNITS)  TO  AVOID  INFINITE  WTS 
DISMIN=Q . 15 
DO  400  1=1,  NCOL 
DO  400  J=l,  NROW 
SUMWT=0 . 0 

DO  100  IT= 1 ,  NUMTOT 

DIST= ( FLOAT( I ) -XG (IT) ) *  *  2+ ( FLOAT( J ) -YG ( IT) )**2 
DIST=SQRT( DIST ) 

IF  (DIST. LE.DISMIN)  DIST=DISMIN 
WT( IT ) = 1 . 0/( DIST* DIST ) 

SUMWT=SUMWT+WT( IT) 

100  CONTINUE 

NORMALIZE  WEIGHTS 

DO  120  IT=1,  NUMTOT 

WT  ( IT)=WT( IT)/SUMWT 
120  CONTINUE 

IF  (J  .Eg.  11)  THEN 

1F( DBUG( 4 ) )  PRINT  9030,  I,  J 

IF( DBUG ( 4 ) )  PRINT  9010,  ( WT( IT ), IT= 1 , NUMTOT ) 

END  IF 

:  MAKE  GRID  POINT  ANALYSIS  USING  WTS  AND  STA  DATA 
DO  350  K=1 ,  NLVL 
U( I , J , K ) =0 . 0 
V( I , J , K ) -0 . 0 


/  „• 
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DO  300  IT=1,  NUMTOT 

U(I,J,K)=U(I,J,K) +WT (IT)*USIG(IT,K) 

V ( I , J , K ) =V ( I , J , K } +WT( IT)*VSIG(IT,K) 

300  CONTINUE 

350  CONTINUE 

400  CONTINUE 

IF  ( DBUG ( 4 ) )  THEN 
DO  430  K=2 ,  NLVL 
PRINT  9035 ,  K 
DO  420  JR=1 ,  NROW 
JP=NROW+l-JR 

420  PRINT  9031,  ( U( IP , JP , K ) , IP= II , 12 ) 

PRINT  9038,  K 
DO  425  JR=1 ,  NROW 
JP=NROW+l-JR 

425  PRINT  9031,  ( V ( IP , JP , K ) , IP= II , 1 2 ) 

430  CONTINUE 

PRINT  9002 
END  IF 

9001  FORMAT  (/  '  BEGIN  SUBROUTINE  GPAN  • /) 

0002  FORMAT  (//'  END  OF  SUBROUTINE  Gl-ni.  , 

9 U 1 0  FORMAT  ( 3X , 8F6 . 1 ) 

9030  FORMAT  (/'  WTS  FOR  STATIONS  FOR  GRID  POINT  X,Y='2I3) 

9031  FORMAT  (1X,21F5.1) 

9035  FORMAT  (/'  U  COMPONENT  AT  LEVEL  = ' 13/ ) 

9038  FORMAT  ('  V  COMPONENT  AT  LEVEL  ='I3/) 

RETURN 

END 

C 


SUBROUTINE  LEVWND( DBUG, IG , ZFLAT) 

c 

C  LINEARLY  INTERPOLATES  WINDS  FROM  SIGMA  SURFACES  TO  HORIZONTAL  PLANES  OR 
C  TO  A  HEIGHT  OF  10  M  ABOVE  THE  LOCAL  TERRAIN  FOR  THE  1ST  LEVEL 
C  — F. LUDWIG,  12/23/85 

C 

LOGICAL  DBUG (15) 

DIMENSION  ZFLAT(6) 

COMMON  /BLHT/  BLT( 25 , 25 ) , HSITE,  AVTHK,  SLFAC, STHK, BLGRX, BLGRY 
COMMON/RARS/RHS ( 25 , 25 , 6 ) 

COMMON/CSFC/SFCHT( 25,25) , SIGMA ( 6 ) 

COMMON/UARS/U ( 25 , 25 , 6 ) , UA( 25 , 25 , 6 ) , V ( 25 , 25 , 6 ) , VA( 25 , 25 , 6 ) 

COMMON /WAR  S/W (  25,25,6), WA (  25, 25,6) 

COMMON  /LIMITS/NCOL , NROW , NLVL, NCOLM1 , NROWM1 , 

$  LOWIX( 5,3), LOWIY (5,3), SFCMAX 

r 

C  STATEMENT  FUNCTION  FOR  LINEAR  INTERPOLATION. 

c 

VALINT(A1 ,A2 , X) =A1+ (A2-A1 ) *X 

IF  ( DBUG (15))  PRINT*,  'START  LEVWND' 

C 

C  GET  LOWEST  PTS  AT  SFC  &  TOP  OF  BOUNDARY  LAYER,  THEN  DEFINE  LAYER  HTS . 

C 

SFCMIN=SFCHT( LOWIX( 1 , IG ) , LOWIY ( 1 , IG ) ) 

THICK=RHS(LOWIX(l, IG) , LOWIY (1, IG) , NLVL) 

DO  25  L=1 , NLVL 

ZFLAT ( L) =SFCMIN+THICK*SIGMA ( L) 

IF  ( DBUG ( 14 ) )  PRINT  *,  ZFLAT( L ) , ’ HT  AT  LEV ' , L 
25  CONTINUE 

DO  100  IX=1 , NCOL 
DO  100  I Y=1 , NROW 

.*  1ST  LEVEL  IS  NOT  FLAT;  IT  IS  SET  AT  10  M  ABOVE  THE  LOCAL  TERRAIN. 

C 

ZFLAT ( 1)=SFCHT( IX, IY)+10.0 


DO  75  IFLAT=1 , NLVL 

IS  THIS  LEVEL  MORE  THAN  1/2  METER  ABOVE  GROUND? 

IF( ZFLAT( IFLAT) -SFCHT{ IX, IY ) .GT.0.5)  THEN 
DO  50  IZ=2 , NLVL 

FIND  THE  FLOW  SFCS  BETWEEN  WHICH  THIS  FLAT  SFC  LIES  5.  DO  LOG  INTERPOL 


IF(RHS( IX, IY, IZ)+SFCHT( IX, IY) .GE. ZFLAT( IFLAT) ) 

$  THEN 

ZBLO=RHS( IX, IY,  12-1) 

IF  (ZBLO.LT.l. )ZBLO=1.0 
DZT=ALOG10 ( RHS  ( IX, IY, IZ)/ZBLO) 

ZFLT=ZFLAT( IFLAT) -SFCHT( IX, IY) 

IF  ( ZFLT  .LT.  1.0)  ZFLT=1.0 
DZF=ALOGl  0  (  ZFLT/ZJiLi  > ) 

RATIO=DZF/DZT 

IF( RATIO. GT. 1 . 001 ) STOP  'BAD  RATIO  IN  LEVWND ' 

UA( IX, IY, IFLAT) =VALINT(U( IX, IY, IZ-1 ) ,U( IX, IY, IZ ) , 

$  RATIO) 

VA( IX, IY, IFLAT )=VALINT(V( IX, I Y , IZ-1 ) , V ( IX , IY , IZ ) , 

$  RATIO) 

WA( IX, I Y , IFLAT ) =VALINT( W ( IX, IY , IZ-1 ) , W( IX, IY , IZ ) , 

$  RATIO) 

GO  TO  75 
END  IF 

>  CONTINUE 
ELSE 

UA( IX, I Y , IFLAT) =0 . 0 
VA( IX, I Y , IFLAT ) =0 . 0 
WA( IX, I Y , IFLAT ) =0 . 0 
END  IF 

>  CONTINUE 
jU  CONTINUE 

IF( DBUG ( 14 ) )  PRINT  *,  'END  LEVWND' 

RETURN 
END 

*********************************************************************** 


SUBROUTINE  REINT ( DBUG , L ) 

INTERPOLATES  LOG-L I NEARLY  TO  GET  U&V  FOR  LEVEL  L  USING  NEXT  LOWER  LEVEL 
(.  THE  TOP  LEVEL.  LOWEST  (.  HIGHEST  LAYERS  ARE  DETERMINED  1ST,  THEN  3,4  ETC 

LOGICAL  DBUG (15) 

COMMON  /LIMITS/NCOL, NROW, NLVL, NCOLM1 , NROWM1 , 

$  LOWIX(  5,3),  LOWI Y  (5,3),  SFCjMAX 

COMMON  /RARS/RHS (25,25,6) 

COMMON  /UARS/U( 25, 25, 6) ,UA( 25,25,6),V(25,25,6),VA(25,25,6) 

STATEMENT  FUNCTION  FOR  INTERPOLATING 

VALU ( XI , X2 , UT, D1 ) =X1+  ( X2-X1 ) *D1/DT 
ZERO=0 . 0 

IF  (DBUG (13))  PRINT* ,' START  REINT' 

DO  100  1=1 , NCOL 
DO  100  J=1 , NROW 

IF( RHS ( I , J , L-l ) . GT .0.0)  THEN 
ZB=RHS( I, J,L-1) 

IF( ZB . LT . 10 . ) ZB= 10 . 0 
DLT=ALOG10 ( RHS ( I , J , NLVL )/ZB ) 

DL1=ALOG10 ( RHS ( I,J,L)/ZB) 

U( I , J , L)=VALU(U( I , J, L-l ) ,U( I , J , NLVL ) , DLT, DL1 ) 

V ( 1 , J , L ) =VALU ( V ( I , J , L-l ) , V ( I , J , NLVL ) . DLT, DL1 ) 


ELSE 

IF( RH3{ X , J , L) . GT .0.5)  THEN 

C 

C  ASSUME  WIND=0  BELOW  0.5  METER 

C 

DLT=ALOG10(RHS( I , J , NLVL)/0 . 5 ) 

DL1=ALOG10 (RHS( I , J , L)/0 . 5 ) 

U( I ,  J, L)=VALU{ ZERO, U( I , J,NLVL) , DLT, DL1 ) 

V ( I , J , L ) =VALU ( ZERO , V ( I , J , NLVL ) , DLT , DL1 ) 

ELSE 

U( I , J, L)=0 . 0 
V( I , J , L ) =0 . 0 
END  IF 
END  IF 

100  CONTINUE 

IF( DBUG ( 13 ) )  PRINT* , ' END  REINT ' 

RETURN 

END 

C 

£*«********»************************************************************** 

SUBROUTINE  RES IG ( DBUG , IG ) 

C 

C  REDEFINES  THE  HEIGHTS  OF  THE  SIGMA  SURFACES  BASED  ON  WIND  SPEED  OVER 
C  THE  LOWEST  TERRRAIN  HEIGHTS  &  LAPSE  RATES  AT  THE  VARIOUS  SIGMA  LEVELS 
C  --THE  UNDERLYING  CONCEPT  IS  SIMILAR  TO  THAT  OF  THE  "CRITICAL  STREAMLINE. 

C 

LOGICAL  DBUG (15) 

DIMENSION  RHSLO( 6 ) , DZMAX( 6 ) , COMPDZ ( 5 ) 

COMMON  /SOUND/  PTLAPS ( 6 ) , TO ( 6 ) 

COMMON  /LIMITS/NCOL , NROW , NLVL , NCOLM1 , NROWM1 , 

$  LOWIX( 5,3), LOWIY (5,3), SFCMAX 

COMMON  /RARS/RHS (25,25,6) 

COMMON  /BLHT/  BLT( 25 , 25 ) , HSITE ,  AVTHK,  SLFAC , STHK , BLGRX , BLGRY 
COMMON  /CSFC/  SFCHT( 25 , 25 ) , SIGMA( 6 ) 

COMMON  /UARS/U{ 25,25,6), UA (25,25, 6),V(25, 25, 6) , VA( 25,25,6) 

DATA  RHSLO, CMPRES  /6*0.0,0.5/ 

C  GETTING  RMS  AVERAGE  WIND  SPEED  OVER  LOW  TERRAIN  GRIDS  &  LOWEST  PT. 

IF  ( DEUG ( 5 } )  PRINT  *(  'START  RESIG' 

LX=LOWIX( 1 , IG ) 

LY=LOWIY ( 1 , IG ) 

SFCMIN=SFCHT( LX, LY ) 

HIRISE=SFCMAX-SFCMIN 
DO  100  1=1,6 
RMSSPD=0. 0 
DO  50  J=l, 5 

JX=LOWIX(J, IG) 

JY=LOWIY ( J, IG) 

RMSSPD=RMSSPD+U( JX, JY, I ) ** 2+V( JX , JY , I ) **2 
50  CONTINUE 

RHSLO ( I ) =RHS ( LX , LY , I) 

SPD=SQRT( RMSSPD/5 .0) 

IF  { DBUG ( 5 ) )  PRINT  *,'RMS  SPEED  AT  LEV',I,SPD 
DTHETA=PTLAPS( I) 

IF  (DTHETA . LT . 1 . E- 20 )  DTHETA= 1 . E- 20 
DZMAX( I )=SPD/SQRT( 9 . 8  * DTH ETA/TO ( I ) ) 

IF( DBUG ( 5 ) )  PRINT* , 'DZMAX= ’ ,DZMAX( I ) , 'AT  LEV ' , I 
100  CONTINUE 

C  FLOW  AMPLITUDES  ARE  LIMITED  TO  TERRAIN  AMPLITUDE  ON  FIRST  PASS 
C 

IF  (DZMAX(NLVL) .GT.HIRISE)  DZMAX ( NLVL ) =HIR ISE 
C  GETTING  1/2  THE  SIGMA  SFC  SEPARATIONS  OVER  THE  LOW  SFC  POINT. 
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DO  125  L=1 , NLVL- 1 

COMPDZ ( L) =CMPRES* ( RHSLO( L+l ) -RHSLO( L) ) 

125  CONTINUE 

C 

C  DO  NOT  LET  SEPARATION  BETWEEN  SFCS  COMPRESS  BY  MORE  THAN  CMPRES  FACTOR 

C 

DO  150  L=NLVL-1 ,1,-1 

IF  ( DZMAX (  L)  .  GT.  ( DZMAX(  I.+  l ) +OOMPDZ  ( I.)  )  )  DZMAX { I. )  =  DZMAX (  L+ 1  ) 

$  f COMPUZ ( L ) 

IF  (DZMAX(L) .GT.  HIRISE)  DZMAX( L ) =HIRISE 
IF(DBUG( 5 ) )  PRINT*, ’DZMAX=' ,DZMAX(L) , ’AT  LEV',L 
150  CONTINUE 

C 

C  NOW,  GO  BACK  &  MAKE  SURE  THAT  HIGHER  SURFACES  DO  NOT  RISE  MORE  THAN 
C  SURFACE  DIRECTLY  BELOW,  IF  THE  SFC  BELOW  IS  ABOVE  THE  TERRAIN  HEIGHT 
C  --SIMILAR  TO  THE  CONSTRAINT  THAT  THE  GREATEST  ALLOWABLE  DEVIATION  IS 
C  TO  PARALLEL  THE  TERRAIN. 

C 

DO  165  L=3 , NLVL 

IF  (  ( RHSLO{ L-l ) +DZMAX( L-l ) }  .LT.  SFCMAX  )  GO  TO  165 
IF  (DZMAX(L) .GT. DZMAX (L-l ) )  DZMAX( L)=DZMAX( L-l ) 

165  CONTINUE 

DO  200  IX=1 , NCOL 
DO  200  IY  =  1 , NROW 
C 

C  LOCAL  RISE  OF  SIGMA  SFC  WILL  BE  PROPORTIONAL  TO  THE  MAXIMUM  RISE  (WHICH 
C  OCCURS  OVER  HIGHEST  TERRAIN  POINT)  AND  THE  RATIO  OF  THE  RISE  OF  THE 
C  LOCAL  GROUND  SURFACE  TO  THE  MAXIMUM  TERRAIN  ELEVATION  DIFFERENCE. 

C 

HERE=SFCHT( IX, IY ) 

ZRATIO= ( HERE-SFCMIN ) /HIRISE 
DO  175  L=1 , NLVL 

RISE=ZRATIO*DZMAX( L) 

RHS( IX, IY, L)=RHSLO( L)+SFCMIN+RISE-HERE 
175  CONTINUE 

C 

C  DEFINE  LOCAL  BOUNDARY  LAYER  THICKNESS 

(7 

BLT(  IX,  I Y )  =RHS(  IX,  IY ,  NLVL)  +  HERE 
200  CONTINUE 

IF  ( DBUG ( 5 ) )  THEN 
DO  300  L=1 , NLVL 
PRINT  6000, L 
DO  250  IY=NROW, 1,-1 

PRINT  6001,  IY , (RHS( IX, I Y , L) , IX=1 , NCOL) 

WRITE  (15,6003)  ( RHS( IX, I Y , L) +SFCHT( IX, IY ) , IX=1 , NCOL) 

250  CONTINUE 

PRINT  6002, (IX, IX=1, NCOL) 

300  CONTINUE 

END  IF 

tOOo  FORMAT ( 1H1 , ' LEVEL= ',12) 

o001  FORMAT( 1H0 ,12, 25F5 . 0 ) 

cOG2  FORMAT ( IX, 'COL' , IX, 25 ( 12 , 3X) ) 

6003  FORMAT(1X,25F5.0) 

IF  ( DBUG ( 5 ) )  PRINT  * , ' END  RESIG' 

RETURN 
END 

*****t*t***tttt*********t***********t*****t**t**********«t»************* 

SUBROUTINE  SETBLT(DBUG) 

’**  THIS  SUBROUTINE  SETS  THE  HEIGHT  OF  THE  BNDY  LAYER  TOP.  AVTHK  IS 
J  AVER.  BL  THICKNESS  OVER  AREA.  SLFAC  CONTROLS  THE  SLOPE,  IF  0  THE 
:  TOP  IS  FLAT,  IF  1  THE  BL  TOP  FOLLOWS  THE  TERRAIN.  HSITE  IS  HT  OF 
C  THE  ANCHOR  POINT  (SITE),  STHK  IS  THE  SMALLEST  BL  THICK  ALLOWED. 


C  BLGRX  IS  B  LYR  HT  GRADIENT  TO  E. 

C  BY  R.  M.  ENDLICH  SRI  INTNL 
C  LAST  REVISION  JUNE  '84 
LOGICAL  DBUG (15) 

DIMENSION  B(25,25),IB(25,25) 

COMMON  /LIMITS/MCOL,  NKOW ,  NLVL,  NCoLMJ. ,  , 

$  LOWIX( 5,3), LOWIY (5,3), SFCMAX 

COMMON/CVOS/  RCM, RMF, IV,DSCRS, IXCRS , JYCRS , IXMED, JYMED, 

$  IXFIN , JYFIN 

COMMON  /SITE/  IXS,  JYS,  THSITE,  IGRID 

COMMON  /BLHT/  BLT{ 25 , 25 ) , HSITE,  AVTHK,  SLFAC , STHK, BLGRX, BLGRY 
COMMON/CSFC/  SFCHT( 25 , 25 ) ,  SIGMA(6) 

COMMON  /PRLIM/  II,  12 
I F  { DBUG ( 6 ) )  PRINT  9001 
THK=AVTHK 

IF  (IGRID. GT.l)  THK=THSITE 
IF  (IGRID  .NE.l)  GO  TO  5 
BLX=BLGRX 
BLY= BLGRY 
IX= IXCRS 
JY= JYCRS 

5  CONTINUE 

IF  (IGRID  .NE.  2)  GO  TO  6 
BLX=  BLGRX/ (RCM) 

BLY=BLGRY/(RCM) 

IX=IXMED 
JY= JYMED 

6  CONTINUE 

IF  (IGRID  . NE . 3 )  GO  TO  7 
BLX=BLGRX/ ( RCM*RMF ) 

BLY=BLGRY/ ( RCM*RMF ) 

IX=IXFIN 

JY=JYFIN 

7  CONTINUE 
ITER=0 

10  ITER=ITER+1 
SUM1=0 . 0 
01=0-0 

DO  50  1=1 , NCOL 

DO  50  J= 1 , NROW 

BLT( I , J ) =THK+ ( SLFAC *SFCHT( I , J ) ) + ( 1 . 0-SLFAC ) *HSITE 
ADD  EASTWARD  GRADIENT  TO  BLT  FROM  COL.  OF  ANCHOR  POINT 
BLT ( I , J ) =BLT ( I ,  J  )  +  (I~IX)*( BLX/FLOATJ ( NCOL ) ) 

C  ADD  NORTHWARD  GRADIENT  FROM  ROW  OF  ANCHOR  POINT 

BLT( I , J ) =BLT( I ,  J )  + ( J-JY ) * (BLY/FLOATJ ( NROW) ) 

IF  ( SFCHT ( I , J ) . GT . ( BLT ( I , J ) -  STHK))  BLT( I , J ) =SFCHT( I , J ) 

+  +STHK 

SUM1=SUM1+(BLT( I,J)-SFCHT( I, J) ) 

Q1=Q1+1 . 0 

50  CONTINUE 

IF( IGRID. GT. 1 )  GO  TO  51 
ATH=SUM1/Q1 
THK=THK+ ( AVTHK -ATH ) 

THSITE=BLT( IXS , JYS ) -SFCHT( IXS, JYS) 

IF( DBUG( 6 ) )  PRINT  9010,  AVTHK, ATH,  THSITE 
DIFF=ABS(AVTHK-ATH) 

IF  ( ITER .GT.  9)  GO  TO  52 

IF  (IGRID  , EQ.  1  .AND.  DIFF.GT.1.0)  GO  TO  10 

51  CONTINUE 

IF( IGRID . EQ . 1 )  GO  TO  52 

v’  MAKE  BL  THICKNESS  AT  ANCHOR  POINT  (THSITE)  THE  SAME  FOR 
MED  AND  FINE  GRIDS  AS  IT  WAS  FOR  COARSE  GRID 
ATH=SUM1/Q1 

THSITE2=BLT{ IXS , JYS ) -SFCHT( IXS, JYS) 

THK=THK+ ( THSITE-THSITE2 ) 

DIFF2=ABS ( THSITE-THSITE2 ) 


IF( DBUG ( 6 } )  PRINT  9010,  AVTHK , ATH , THSITE2 
IF( ITER . GT . 9 )  GO  TO  92 
IF(DIFF2.GT. 1.0)  GO  TO  10 
52  CONTINUE 

DO  55  JP=1 , NROW 
DO  55  IP=1 , NCOL 

B(  IP, JP)  =  BLT( IP, JP) 

IB( IP , JP ) =JNINT( B( IP , JP ) ) 

55  CONTINUE 

IF( DBUG ( 6 ) )  PRINT  9115,IGRID 
DO  60  JP=1 , NROW 
JR=NROW+  1  -JP 

60  IF(DBUG(6) )  PRINT  9105,  { IB( IP , JR ) , IP=I1 , 12 ) 

DO  65  JP=l,NROW 
DO  65  1P=1 , NCOL 

B( IP, JP)=(BLT( IP, JP)-SFCHT(IP, JP) ) 

IB( IP , JP ) =JNINT( B( IP , JP ) ) 

65  CONTINUE 

I F ( DBUG ( 6 ) )  PRINT  9020 
DO  70  JP=1 , NROW 
JR = NROW  "tl-JP 

70  IF(DBUG(6) )  PRINT  9105,  ( IB( IP, JR ) , IP=I1, 12 ) 

IF(DBUG( 6 ) )  PRINT  9002 

9001  FORMAT  ( 1H1 , '  BEGIN  SUBROUTINE  SETBLT'/) 

9002  FORMAT  (//•  END  OF  SUBROUTINE  SETBLT’/) 

9010  FORMAT  (’  INITIAL  AV.  THICKNESS,  M=’F10.1, 

+  ’  ACTUAL  AV.  THICKNESS  =’F10.1, 

+  '  SITE  THICKNESS  ='F10.1/) 

9020  FORMAT ( 1H1, ’ BNDY  LAYER  THICKNESS  IN  M’/) 

9105  FORMAT  (/IX, 2515) 

9115  FORMAT  (1H1,’  HEIGHT  OF  BNDY  LAYER  TOP,  M,  GRID  =  ’I3) 
RETURN 
END 

it*********************************************************** 


SUBROUTINE  SETMAT ( VALUE , ARRAY , NUM1 , NUM2 ) 

INITIALIZES  ALL  ELEMENTS  OF  ARRAY  TO  VALUE 
REVISON  SEPT.  1978 

DIMENSION  ARRAY ( NUM1 , NUM2 ) 

DO  10  1=1, NUM1 
DO  10  J=1 , NUM2 

ARRAY ( I , J ) = VALUE 
CONTINUE 
RETURN 
END 

£»*****»*********»***********************«******«* 

SUBROUTINE  STRAT( DBUG , IG ) 

V. 

0  READS  SOUNDING  INFORMATION , CALCULATES  LAPSE  RATES  &  OTHER  PARAM 
C  ETERS  FOR  VARIOUS  SIGMA  LEVELS 

C 

LOGICAL  DBUG( 15 ) 

DIMENSION  Z( 50) ,T( 50) ,P( 50) , DPTDZ(49) , ZMID(49) 

COMMON  /RARS/RHS (25,25,6) 

COMMON  /LIMITS/NCOL, NROW, NLVL, NCOLM1 , NROWM1 , 

$  LOWIX( 5,3), LOWIY (5,3), SFCMAX 

COMMON  /SOUND/  PTLAPS ( 6 ) , TO( 6 ) 

C  STATEMENT  FUNCTION  FOR  POTENTIAL  TEMPERATURE 


^  *  *  *  *  » 


10 


THETA (P,T)=(T+273 . 13) * ( ( 1000 . /P )•* 0 . 288 ) 


READ  NO.  OF  HEIGHTS  IN  SOUNDING  &  DATA  INPUT  FORM: 

ITYP=1-- 

HEIGHT, POT  TEMP  LAPSE  RATE 

ITYP=2  — 

HEIGHT,  TEMP  &  PRESSURE 

READ( 13 , 6001 )  NHITES, ITYP 
READ(13,6002) (Z(I) ,T{ I ) ,P( I) , 1=1, NHITES) 

IF  ( DBUG ( 8 ) )  PRINT  *, 'START  STRAT ' 

READ( 13 , * )  NHITES, ITYP 

IF  ( DBUG { 8 ) )  PRINT  *,’ NHITES, ITYP ' ,  NHITES, ITYP 
DO  10  1=1, NHITES 

READ( 13 , * )  Z{ I) ,T( I) , P( I) 

IF  { DBUG ( 8 ) )  PRINT  *,  Z(I ) ,T( I) ,P( I) 

10  CONTINUE 

IF(ITYP.EQ.l)  THEN 
DO  20  1=1 , NHITES-1 
DPTDZ ( I ) =P ( I ) 

ZMID( I ) =  Z ( I ) 

20  CONTINUE 

ELSE 

DO  30  1=1, NHITES-1 

DPTDZ ( I ) = ( THETA ( P ( I + 1 ) ,T(I+1) )-THETA(P(I) ,T(I) ) )/ 

S  (Z(I+1)-Z(I)) 

ZMID( I)=0.5*(Z(I)+Z(I+1)J 
30  CONTINUE 

END  IF 

2  GETTING  AVERAGE  HEIGHT  OF  LOWEST  SIGMA  SFCS 

DO  100  1=1,6 
ZBAR=0 . 0 

2  GET  AVERAGE  HT.  ABOVE  LOW  SPOTS  FOR  SIG  SURFACES 
DO  40  J=1 , 5 

ZBAR=ZBAR+RHS( LOWIX( J , IG ) , LOWIY ( J , IG ) , I ) 

40  CONTINUE 

ZBAR=ZBAR/5 . 0 
DO  50  J=2, NHITES 

IF ( ZMID( J ) .GE.ZBAR. OR. J.EQ. NHITES)  THEN 
PTLAPS( I ) =DPTDZ ( J-l ) + { ZBAR-ZMID( J-l ) ) * 

$  ( DPTDZ (J ) -DPTDZ (J-l ) )/( ZMID( J )-ZMID( J-l) ) 

TO( I )=T(J-1)+(ZBAR-ZMID( J-l) )* 

$  (T(J)-T(J-l) }/(ZMID(J)-ZMID(J-l) ) 

TO( I ) =TO( I )  +  273 . 13 
GO  TO  75 
END  IF 

5u  CONTINUE 

-'5  IF  ( DBUG  (  8  )  )  PRINT*,  1  PTLAPS  AT  LEV’  ,  I,  ’  =  ’  ,PTLAPS(  I  ) 

100  CONTINUE 

IF  ( DBUG ( 6 ) )  PRINT*, ’END  STRAT’ 
tOOl  FORMAT  (213) 

0002  FORMAT  (3F10.3) 

RETURN 

END 

C 

Citititttittiiututiintttiotitioitttitmttmtiitiittiitxtttttiii 


SUBROUTINE  TOPO ( NUM , DBUG ) 

2  READ  AND  COMPUTE  TOPOGRAPHY  AT  GRID  POINTS 
C  LAST  REVISION  OCTOBER  ’84. 

C  IF  NUM=0  READS  TERRAIN  HEIGHTS  FOR  ALL  GRIDS. 
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onn  or ;  r. 


IF  NUM  .GT.  0  PICKS  TERRAIN  HTS  ,  CALLS  SETBLT  TO  ESTABLISH 
BNDY  LYR  TOP,  AND  COMPUTES  RELATIVE  HTS  ( RHS )  FOR  PROPER  GRID. 

LOGICAL  DBUG (15) 

COMMON  /LIMITS/NCOL, NROW, NLVL, NCOLM1 , NROWM1 , 

$  LOWIX( 5,3), LOWIY (5,3), SFCMAX 

COMMON /CTOP/  MCRS , NCRS , MMED , NMED , MFIN , NFIN , NGR ID 
COMMON/RARS/RHS( 25,25,6} 

COMMON/CSFC/SFCHT( 25,25) , SIGMA ( 6 ) 

COMMON  /BLHT/  BLT( 25 , 25 ) , HSITE,  AVTHK,  SLFAC, STHK , BLGRX, BLGRY 
COMMON  /PRLIM/  II,  12 

DIMENSION  HTCRS( 25,25),  HTMED( 25 , 25 ) ,  HTFIN(25,25) 

TO  ACCOUNT  FOR  STABLE  FLOWS,  THE  LOWEST  SIGMA  SFCS  SHOULD 
INTERSECT  HIGH  TERRAIN.  THE  LIMIT  FOR  THIS  INTERSECTION  IS 
TERLIM  (IN  M) . 

IF(DBUG( 9 } )  PRINT  9001, NUM 
IF  (NUM.GT. 0)  GO  TO  10 
C  READ  TERRAIN  HEIGHT  VALUES  AT  GRID  POINTS  IN  METERS,  ALL  GRIDS 
IF ( DBUG ( 9 ) )  PRINT  9006 
IF( DBUG{ 9 ) )  PRINT  9003 

C  IN  HEIGHT  DATA  FILE  NORTHERN  ROW  IS  FIRST,  SO  INVERT  ORDER. 

C  READ  AND  PRINT  HEIGHTS  AT  COARSE  GRID  POINTS 
DO  8  J=1 , NCRS 
JR=NCRS+1-J 

READ( 11 , * )  ( HTCRS( I , JR ) , 1=1 , MCRS ) 

8  CONTINUE 

IG=1 

IF  ( DBUG ( 9 ) )  THEN 
DO  118  J=1 , NCRS 
JR=NCRS+1-J 

PRINT  4,  ( HTCRS( I , JR ) , 1=11, 12) 

118  CONTINUE 
END  IF 

C  READ  AND  PRINT  MEDIUM  GRID  HEIGHTS 
IF  (NGRID.LT. 2)  GO  TO  120 
IF( DBUG ( 9 ) )  PRINT  9004 
DO  9  J=1 , NMED 
JR=NMED+ 1- J 

READ( 11 , * )  ( HTM£D( I , JR ) , 1=1 , MMED ) 

9  CONTINUE 

IG=2 

IF  ( DBUG ( 9 ) )  THEN 
DO  119  J=1 , NMED 
JR=NMED+1-J 

PRINT  4,  (HTMED( I , JR) , 1=1, MMED) 

119  CONTINUE 
END  IF 

120  CONTINUE 

C  READ  AND  PRINT  FINE  GRID  HEIGHTS 
IF  (NGR ID  .NE.  3 )  GO  TO  214 
DO  210  J=1 , NFIN 
JR=NFIN  +1-J 

READ( 11 , * )  ( HTFIN ( I , JR ) , 1=1 , MFIN ) 

210  CONTINUE 

I F ( DBUG ( 9 ) )  PRINT  9005 
IG=3 

IF  ( DBUG ( 9 ) )  THEN 
DO  212  J= 1 , NFIN 
JR=NFIN+ 1- J 

PRINT  4,  ( HTFIN ( I , JR } , 1=1 , MFIN ) 
zU  CONTINUE 

END  IF 

/ 1 4  CONTINUE 
GO  TO  150 

10  CONTINUE 

DO  15  J= 1 , NROW 


DO  15  1=1 ,NCOL 

IF  (NUM  .NE.  1)  GO  TO  11 

SFCHT( 1 , J ) =HTCRS( I , J ) 

11  CONTINUE 

IF  (NUM  .NE.  2)  GO  TO  12 
SFCHT( 1 , J } =HTM£D( I , J ) 

12  CONTINUE 

IF  (NUM  . NE. 3 )  GO  TO  13 
SFCHT( I ,  J  ) =HTFIN ( I ,  J ) 

13  CONTINUE 
15  CONTINUE 

C  FIND  HIGHEST  &  LOWEST  TERRAIN  POINTS 

C 

CALLXTR EM ( N UM , DBUG ) 

C**  SET  BNDY  LAYER  TOP  (ARRAY  BLT) . 

CALL  SETBLT(DBUG) 

C  DENOTE  GEOMETRIC  HEIGHT  ABOVE  TERRAIN  BY  RHS 
DO  67  J=1 , NROW 
DO  67  1=1, NCOL 

ZVAR=BLT(I,J)-SFCHT(I,J) 

DO  67  K=1 , NLVL 

RHS ( I , J, K ) =SIGMA( K ) *ZVAR 
67  CONTINUE 

150  IF( DBUG ( 9 ) )  PRINT  9002 
2  FORMAT( 8F10 . 2 ) 

4  FORMAT(/ , 5X, 21F5 . 0 ) 

3001  FORMAT  (/'  BEGIN  SUBROUTINE  TOPO,  NUM='I3) 

9002  FORMAT  (/’  END  OF  SUBROUTINE  TOPO’/) 

^003  FORMAT  (1H1,’  TERRAIN  HTS,  COARSE  GRID,  METERS1/) 

9004  FORMAT  (1H1,1  TERRAIN  HTS,  MEDIUM  GRID,  METERS’/) 

3005  FORMAT  (1H1,‘  TERRAIN  HTS,  FINE  GRID,  METERS’) 

9006  FORMAT  (/’  PRINTOUT  IS  REVERSE  OF  INPUT  -  HAS  NORTH  ROW  1ST’/) 
RETURN 
END 

C 

(j*************t**********a********************t**«*********************** 

C 

SUBROUTINE  WXANAL( NUM, DBUG ) 

C  THIS  SUBROUTINE  DOES  THE  FOLLOWING.  WHEN  NUM=IGRID=1  IT 
C  READS  IN  DATA  FROM  WIND  PROFILES  AND  ASSIGNS  WINDS  TO 
C  SIGMA  LEVELS,  READS  NWS  STATIONS,  COMPUTES  THE  GEOSTROPHIC 
C  WIND  AND  ASSIGNS  WINDS  TO  SIGMA  LEVELS,  COMBINES  SOUNDINGS  AND 
C  HOURLY  DATA,  AND  MAKES  INITIAL  OBJECTIVE  WIND  ANALYSIS  USING 
C  A  WT  FACTOR  INVERSELY  PROPORTIONAL  TO  DISTANCE  SQUARED. 

WHEN  NUM  .GT.  1  (IGRID=1  OR  2 )  IT  SELECTS  INITIAL  WINDS  FROM 
C  THE  NEXT  COARSER  GRID. 

C  BY  RM  ENDLICH,SRI  INTN’L,  JAN  ’85 


LOGICAL  DBUG (15), KEY 

COMMON  /LIMITS/NCOL, NROW, NLVL, NCOLM1 , NROWM1 , 

$  LOWIX( 5,3), LOWIY (5,3) , SFCMAX 

COMMON  /CVOS/  RCM,KMF, IV,DSCRS, IXCRS, JYCRS, IXMED, JYMED, 

+  IXFIN , JYFIN 

COMMON/RARS/RHS ( 25 , 25 , 6 ) 

COMMON  /SITE/  IXS,  JYS,  THSITE ,  I GRID 

COMMON/UARS/U( 25, 25, 6) ,UA( 25, 25, 6) ,V(25, 25,6) , VA( 25, 25, 6) 
COMMON  /NUMOBS/  NUMDOP ,  NUMNWS,  NUMTOT 
DIMENSION  UTEMP(25, 25) ,  VTEMP(25,25) 

I F ( DBUG (10))  PRINT  9001,  NUM 

LINES  TO  STATEMENT  300  ARE  USED  TO  READ  SOUNDINGS  ONLY  FOR  COARSE 
GR ID--NUM= 1 . 

IF  (NUM.GT.l)  GO  TO  200 
KEY= .TRUE. 
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C  READ  WIND  SOUNDINGS  AND  ASSIGN  TO  SIGMA  LEVELS 
CALL  DOPSIG(DBUG,KEY) 

C  READ  NWS  REPORTS  AND  INTERPOLATE  TO  SIGMA  SFCS 
READ{ 1 2 , 9014 )  NUMNWS 
READ( 12 , * )  NUMNWS 
I F ( DBUG (10))  PRINT  9027,  NUMNWS 
IF  ( NUMNWS. GE.l)  CALL  GEOSIG(DBUG) 

C  TOTAL  NUMBER  OF  STATIONS  ( NUMTOT) =NUMDOP+NUMNWS 
NUMTOT=NUMDOP+NUMNWS 
C  MAKE  GRID  POINT  ANALYSIS  OF  DATA 
CALL  GPAN(DBUG) 

C  READ  (12,9014)  NTSOND 

READ  (12,*)  NTSOND 

IF  ( DBUG ( 10 ) )  PRINT  *,  'NO. OP  TSONDE= ', NTSOND 
IF  (NTSOND  . EQ . 1 )  THEN 

IF  ( DBUG (10))  PRINT*,  1 WXANAL  CALL  STRAT, IGR1D= ' , IGR ID 
CALL  STRAT  (DBUG, IGR ID) 

CALL  RESIG(DBUG, IGRID) 

KEY = .FALSE. 

CALL  DOPSIG( DBUG, KEY) 

CALL  GPAN(DBUG) 

END  IF 

DO  50  J=l,  NROW 
DO  50  1=1,  NCOL 
DO  40  LV=2 , NLVL 

C  TO  USE  RHS  NEGATIVE  (BELOW  TERRAIN)  MAKE  WINDS  0. 

IF  ( RHS (I,J,LV) -GE.0.0)  GO  TO  40 
U( I, J,LV)=0. 0 
V  ( I ,  J , LV ) =  0 . 0 
40  CONTINUE 
50  CONTINUE 
GO  TO  300 
200  CONTINUE 

C  SELECT  WINDS  FOR  SMALLER  GRID  FROM  LARGER  GRID 
C  ASSIGN  WINDS  TO  MED.  GRID  FROM  COARSE  GRID 
IF  (NUM  .NE.  2)  GO  TO  216 
DO  215  K=2 , NLVL 
DO  210  1=1, NCOL 
DO  210  J  =  l,NROY< 

IC=IXCRS+JNINT( FLOAT ( I-IXMED)/RCM) 

JC= JYCRS+JN INT ( FLOAT ( J -JYMED ) /R CM ) 

C  FILL  IN  NONZERO  WINDS  FROM  NEAREST  POINTS 


IF 

( U( IC , JC , K ) 

.EQ. 

0. 

.AND.  V ( IC , JC , K ) 

■EQ. 

0.  ) 

IC=IC+1 

IF 

(U(IC,JC,K) 

.EQ. 

0. 

•AND.  V( IC , JC , K ) 

.  EQ. 

0.  ) 

IC=IC-2 

IF 

( U( IC , JC , K ) 

.EQ. 

0. 

.AND.  V ( IC , JC , K ) 

■EQ. 

o. ) 

JC=JC+1 

IF 

( U( IC , JC , K ) 

■  EQ. 

0. 

■AND.  V ( IC , JC, K ) 

.EQ. 

0.  ) 

JC=JC~ 2 

UTEMP( I , J)=U( IC, JC,K) 
VTEMP ( I , J ) =V ( IC , JC , K ) 


<:10  CONTINUE 

DO  212  1=1, NCOL 
DO  212  J= 1 , NROW 
U ( I , J , K ) =UTEMP ( I , J ) 

V( I , J , K ) =VTEMP ( I , J ) 

IF  ( RHS ( 1 , J , K ) . LE .0.0)  U ( I , J , K ) =0 . 0 
IF  ( RHS ( 1 , J , K ) . LE .0.0)  V(I,J,K)=0.0 
212  CONTINUE 

215  CONTINUE 

216  CONTINUE 

C  SELECT  FINE  GRID  WINDS  FROM  MEDIUM  GRID  WINDS 
IF  (NUM  .NE.  3)  GO  TO  300 
DO  230  K=  2 , NLVL 
DO  225  1=1, NCOL 
DO  225  J=1 , NROW 

IC=IXM£D+JNINT( FLOAT ( I-IXFIN )/RMF) 
JC=JYMED+JNINT(FLOAT( J-JYFIN )/RMF) 

'  FILL  IN  NONZERO  WINDS  FROM  NEAREST  POINTS 
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IF 

(U( IC, JC, K) 

•  E Q. 

0. 

.AND.  V ( IC , JC , K ) 

■  EQ. 

o. ) 

IC=IC+1 

IF 

(U( IC, JC,K) 

■  EQ . 

0. 

.AND.  V( IC, JC, K ) 

.  EO  • 

0.) 

IC=IC-2 

IF 

(U(IC,JC,K) 

■EQ. 

0. 

.AND.  V(IC,JC,K) 

.EQ. 

o. ) 

JC=JC+1 

IF 

(U(IC,JC,K) 

•  EQ. 

0. 

.AND.  V(IC,JC,K) 

.EQ. 

o. ) 

JC=JC-2 

UTEMP (I/J)=U(IC,JC,K) 

VTEMP(I,J)=V(IC,JC,K) 

225  CONTINUE 

DO  227  1  =  1 , NCOL 
DO  227  J=1 , NROW 
U  (  I ,  J  ,  K  )  =UTEMP  (  I ,  J  ) 

V  ( I ,  J  ,  K ) = VTEMP ( I ,  J ) 

IF  ( RHS ( I ,  J ,  K ) . LE .0.0)  U(I,J(K)  =  0.0 
IF  ( RHS ( I , J , K ) . LE .0.0)  V(I,J,K)=0.0 
227  CONTINUE 
230  CONTINUE 
3  00  CONTINUE 

I F ( DBUG (10))  PRINT  9002 

^001  FORMAT  (/'  BEGIN  SUBROUTINE  WXANAL,  NUM='I3/) 

9002  FORMAT  (/’  END  OF  SUBROUTINE  WXANAL'/) 

9014  FORMAT  (3X, 15, 2F8 . 2) 

9027  FORMAT  (/'  NUMBER  OF  HOURLY  SURFACE  REPORTS  ='I3/) 
RETURN 


END 


SUBROUTINE  XTREM  ( IG , DBUG ) 

FINDS  5  LOWEST  SFC  ELEVATIONS  &  GRID  PT.  WITH  HIGHEST  ELEVATION 

LOGICAL  DBUG (15) 

DIMENSION  Z LOW ( 5 ) 

COMMON/C SFC/ 5 FCHT{ 25 , 25 ) , SIGMA( 6 ) 

COMMON  /L I MITS/NCOL , NROW , NLVL , NCOLM1 , NROWM1 , 

S  LOWIX( 5,3), LOWI Y (5,3), SFCMAX 

DO  10  1=1,5 

ZLOW ( I ) =99999 . 

10  CONTINUE 

SFCMAX=- 999999. 

DO  100  I Y= 1 , NROW 
DO  100  IX=1 , NCOL 

IF( SFCHT ( IX, IY) . LT. ZLOW( 5 ) )  THEN 
DO  50  1=1,5 

IF( SFCHT ( IX, IY ) . LT . ZLOW( I ) )  THEN 
I F ( I . LT . 5 )  THEN 

DO  45  J=5, If  1, -1 
ZLOW ( J ) =ZLOW( J- 1 ) 

LOWIX( J , IG ) =LOWIX ( J- 1 , IG) 

LOWIY ( J , IG ) =LOWI Y ( J- 1 , IG) 

CONTINUE 
END  IF 

ZLCW( I ) =SFCHT( IX, IY) 

LOW I X ( 1 , 1G)  =  IX 
LOW I Y ( I , IG ) = IY 
GO  TO  75 
END  IF 

50  CONTINUE 

END  IF 

IF(SFCHT( IX, IY) .GT. SFCMAX)  SFCMAX=SFCHT( IX, IY ) 

1 90  CONTINUE 

IF  ( DBUG (11))  PRINT  6000,  IG , ( LOWIX ( I , IG ) , LOWI Y ( I , IG ) , 

$  ZLOW(I) , 1=1,5) 

IF  ( DBUG (11))  PRINT  *,'HIGHEST  ELEVATION= ', SFCMAX 
<  juO  FORMAT  (IX, 'IN  XTREM,  GRID=',I2,’  LOWEST  COL, ROW  HT S’, 

$  5 ( / , 214 , FS . 1  )  ) 

RETURN 
END 
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I  INTRODUCTION 


The  intent  of  this  appendix  is  to  provide  a  review  of  some  recent 
developments  in  the  mathematical  description  of  the  state  of  the  atmos¬ 
phere,  and  of  atmospheric  processes.  The  emphasis  is  on  physical  inter¬ 
pretation  of  the  mathematical  concepts.  Originally,  the  review  was  to 
be  limited  to  fractals  and  fractal  dimension.  The  emphasis  remains 
there,  but  the  concept  of  fractal  dimension  appears  not  only  in  connec¬ 
tion  with  direct  description  of  spatial  and  temporal  distributions  of 
atmospheric  parameters,  but  also  as  a  descriptor  of  equations  used  to 
describe  atmospheric  processes.  It  is  in  connection  with  the  latter 
application  that  the  concept  of  "strange  attractor"  enters  the  picture. 
This  concept  will  be  discussed,  but  not  emphasized. 

To  date,  most  studies  have  focused  on  descriptive  and  interpreta¬ 
tive  applications  of  the  concepts.  Very  little  has  been  done  toward 
incorporating  the  concepts  into  atmospheric  modeling.  This  is  not 
surprising  because  it  is  little  more  than  a  decade  ago  that  the  terms 
"fractal"  and  "fractal  dimension"  were  coined  (Mandelbrot,  1975).  This 
neologism  arose  when  it  was  recognized  that  what  had  appeared  to  be  a 
very  abstract  branch  of  mathematics  could  serve  as  an  analog  for  many 
natural  phenomena.  One  of  the  natural  phenomena  to  which  the  concept 
seemed  most  applicable  was  atmospheric  turbulence.  With  this  relatively 
short  history,  it  is  not  surprising  that  the  concepts  have  yet  to  be 
incorporated  into  fluid  flow  modeling.  I  have  attempted  to  go  somewhat 
beyond  the  usual  bounds  of  a  literature  review  and  have  made  some  tenta¬ 
tive  suggestions  about  how  these  new  concepts  might  be  incorporated  into 
the  flow  modeling  process. 


II  BACKGROUND 
"Unpredictable"  Determinism 
Lorenz  (1963)  observed: 

"certain  hydrodynamical  systems  exhibit  steady-state  flow 
patterns,  while  others  oscillate  in  a  regular  periodic 
fashion.  Still  others  vary  in  an  irregular,  seemingly  hap¬ 
hazard  manner,  and,  even  when  observed  for  long  periods  of 
time  do  not  appear  to  repeat  their  previous  history..." 

"Lack  of  periodicity  is  very  common  in  natural  systems 
and  is  one  of  the  distinguishing  features  of  turbulent  flow. 
Because  instantaneous  turbulent  flow  patterns  are  so  irregu¬ 
lar,  attention  is  often  confined  to  the  statistics  of  turbu¬ 
lence  which,  in  contrast  to  details  of  turbulence,  often 
behave  in  a  regular  well-organized  manner.  The  short-range 
weather  forecast,  however,  is  forced  willy-nilly  to  predict 
the  details  of  the  large  scale  turbulent  eddies--cyclones  and 
anticyclones — which  continually  arrange  themselves  into  new 
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patterns.  Thus,  there  are  occasions  when  more  than  the  sta¬ 
tistics  of  irregular  flow  are  a  very  real  concern." 

Lorenz  (1963)  then  went  on  to  show  that  a  system  could  be  fully 
deterministic  in  the  sense  that  its  state  at  one  step  fully  determined 
its  state  at  the  next  step,  but  would  still  be  essentially  unpredictable 
beyond  a  few  steps,  because  of  extreme  sensitivity  of  the  system  to 
initial  conditions.  Effectively,  these  sensitivities  make  it  impossible 
to  specify  the  state  of  such  a  system  sufficiently  well  that  its  future 
state  can  be  predicted.  The  Navier-Stokes  equations  are  such  a  system. 
By  extension,  those  fluid  systems  such  as  the  atmosphere  that  can  be 
described  by  the  Navier-Stokes  equations  will  be  generally  unpredictable 
in  their  details  and  even  in  their  major  features  beyond  some  limited 
time  in  the  future.  The  deterministic  nature  of  the  equations  and 
pictures  of  turbulent  flows  suggest  that  although  future  states  of  the 
system  are  not  predictable,  they  may  not  be  wholly  without  order.  Thus, 
one  is  encouraged  to  seek  descriptions  of  the  system  that  do  not  obscure 
the  order  and  organization  that  is  present.  The  concepts  of  fractal 
dimension  and  strange  attractors  seem  to  provide  improvements  on  conven¬ 
tional  statistical  measures.  They  augment  the  statistical  measures 
rather  than  replace  them. 


B.  Physical  Space  and  Phase  Space 

In  order  to  illustrate  his  point,  Lorenz  (1963)  used  a  simple  set 
of  equations  to  describe  a  particle  trajectory  governed  by  processes 
closely  related  to  those  causing  thermally  driven  convection.  Those 
equations  are 

X'  =  -oX  +  oY 

Y'  =  -XZ  +  rX  -Y 

Z'  =  XY  -bZ  .  ( II- 1 ) 

Here,  (  ) '  represents  a  rate  of  change  with  respect  to  a  dimensionless 

time;  X,  Y  and  Z  are  locations  and  the  other  parameters  are  constants 

that  can  be  related  to  thermodynamic  properties.  The  above  set  of 

equations  has  come  to  be  called  the  "Lorenz  attractor."  Lorenz  (1963) 

analyzes  the  properties  of  the  equations  to  show  that  certain  choices  of 
constant  values  and  initial  conditions  lead  to  steady-state  or  periodic 
trajectories.  However,  the  most  interesting  cases  are  those  that  are 
neither  stationary  nor  periodic.  He  used  numerical  integration  to  study 
the  subspace  in  which  such  trajectories  are  ultimately  confined. 
Lorenz'  analysis  indicates  that  the  trajectories  become  confined  to  a 
convoluted  surface  that  has  zero  volume.  If  the  trajectories  were 
periodic  (i.e.  they  passed  more  than  once  through  the  same  point,  and 
because  of  the  deterministic  nature  of  the  system,  they  repeated  their 
passage  through  all  subsequent  points),  then  the  solutions  would  be 
confined  to  a  closed,  distorted  curve  that  had  neither  volume  nor  area. 


The  Lorenz  trajectories  will  often  tend  to  cluster  around  a  limited 
number  of  points,  referred  to  as  attractors.  The  trajectory  may  stay 
for  a  number  of  steps  in  the  vicinity  of  one  of  the  attractors  and  then 
make  a  relatively  rapid  leap  into  the  domain  of  another,  stay  there  for 
awhile  and  then  return.  The  nature  of  the  trajectory-containing  sui — 
faces,  and  the  factors  that  govern  shifts  from  one  attractor  to  another 
are  of  considerable  interest  because  of  their  relevance  to  our  under¬ 
standing  of  the  solutions  to  these  systems. 

In  the  case  of  the  Lorenz  attractor,  the  motion  of  the  point 
described  by  the  equations  takes  place  in  three-dimensional  physical 
space.  The  x,  y,  z  coordinates  and  the  u,  v,  w  components  of  motion 
could  be  considered  to  be  descriptors  of  the  "state"  of  the  particle. 
It  is  common  to  represent  the  state  of  a  fluid  system  by  the  values  of 
relevant  parameters  (e.g.  wind  components,  temperature)  at  a  number  of 
grid  points.  Then,  the  evolution  of  the  fluid  flow  can  be  considered  as 
a  trajectory  in  this  high-dimensional  state  (or  phase)  space.  The 
subspace  occupied  by  the  trajectory  will  generally  be  of  lower  topo¬ 
logical  dimension  than  the  state  space  itself.  The  degree  to  which  this 
subspace  fills  the  complete  space  is  determined  by  the  nature  of  the 
equations.  In  many  cases  there  will  preferred  regions  of  the  phase 
space,  i.e.  attractors. 

At  least  intuitively,  the  above  discussion  suggests  that  it  may  be 
possible  to  approximate  fluid  flows  with  fewer  descriptors,  and  that 
descriptions  of  the  shapes  of  the  subspaces  containing  the  solutions 
might  provide  a  basis  for  statistical  descriptions  of  the  flow.  Fractal 
dimension,  as  will  be  discussed  later,  provides  a  measure  of  the  degree 
to  which  the  phase  space  is  filled  by  the  solutions.  It  is  this  fact 
that  has  provided  one  of  the  motivations  for  applying  the  concept  of 
fractal  dimension  to  the  study  of  turbulence.  The  hope  prevails  that 
fractals  will  provide  a  means  for  deriving  statistical  properties  of 
turbulent  fields,  especially  those  associated  with  the  observed  inter- 
mittency  of  most  turbulence. 

There  are  other  reasons,  based  on  physical  arguments,  for  believing 
that  fractal  concepts  show  promise  for  providing  quantitative  descrip¬ 
tions  of  turbulent  fields.  These  other  reasons  will  be  discussed  later, 
but  it  is  worth  noting  here  that  the  physically  descriptive  part  of  the 
theory  can  be  used  directly  to  generate  "realistic"  distributions  of 
parameter  values  in  physical  space,  e.g.  the  shapes  of  surfaces  of 
constant  value  (iso-surfaces)  can  be  generated.  Through  appropriate 
transformations ,  it  is  also  possible  to  estimate  fractal  dimension  from 
a  time  series  of  parameter  values  at  a  single  point. 

This  review  emphasizes  applications  of  fractals  to  distributions  in 
physical  space,  because  the  ultimate  objective  is  to  use  the  concepts  to 
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The  terms  "state  space"  and  "phase  space"  are  both  found  in  the  literature, 
but  the  latter  seems  more  common. 


simulate  small  scale  distributions  from  information  contained  in  larger 
scale  distributions.  In  principle  at  least,  there  is  potential  for 
providing  a  means  of  closure  for  the  large  scale  fluid  flow  equations. 


C.  Statistical  Descriptions 


It  is  natural  to  invoke  statistical  descriptions  when  trying  to 
describe  some  deterministically  unpredictable  behavior.  One  of  the  more 
common  statistical  descriptors  of  turbulent  behavior  is  the  probability 
distribution  associated  with  differences  in  parameter  value  over  a 
specified  distance,  i.e. 


F(q,Ar)  =  Pr[Ac(Ar)  >  q] 


( 1 1— 2 ) 


where  Pr[r  >  s]  is  the  probability  that  the  argument  inequality  is  true; 
Ac(x)  is  the  absolute  value  of  the  difference  in  c  at  two  points  sepa¬ 
rated  by  a  distance  x.  If  the  variable  c  is  uniformly  distributed,  then 
Ac  will  always  be  zero  as  will  Pr[Ac  >  0]. 

For  many  practical  purposes,  it  is  the  tail  of  the  distribution 
that  is  of  greatest  interest.  The  tail  involves  the  large,  infrequent 
fluctuations.  The  nature  of  the  tail  of  the  distribution  will  have 
important  effects  on  the  higher  moments  of  the  probability  function. 
Probability  distributions  will  vary  according  to  the  separation  distance 
between  points.  For  example,  it  is  natural  to  assume  that  the  likeli¬ 
hood  of  exceeding  some  specified  difference  value  increases  as  the 
separation  between  the  points  increases.  The  next  section  describes  a 
simple  functional  relationship  between  the  probability  distributions  for 
different  separations  that  has  been  found  to  be  applicable  to  many 
natural  phenomena. 


D.  Scaling 

Qualitatively,  the  fluctuations  in  a  field  that  is  "scaling"  look 
the  same  regardless  of  magnification.  That  is,  the  large  scale  fluctua¬ 
tions  are  qualitatively  the  same  as  the  middle  scale  fluctuations 
embedded  within  them.  Those  middle  scale  features  are  in  turn  similar 
to  the  still  smaller  fluctuations  superimposed  upon  them.  For  a  scaling 
field  of  the  parameter  c,  this  can  be  quantitatively  defined  as  follows. 
First, 


Ac(Az)  =  c(zQ  +  Az)  -  c(zQ) 


( 1 1-3 ) 


i.e.  Ac  is  the  change  in  c  between  zQ  and  zQ  +  Az.  If  the  field  is 
scaling,  then  the  relationship  between  the  probability  distributions  for 
the  large  scales  and  those  for  the  smaller  scales  is  defined  as  follows: 


Pr[ Ac ( AAz )  >  q]  =  Pr[AnAc(Az)  >  q] 


(II-iJ) 


Scaling  fields  are  fractals 


The  potential  importance  of  fractals  arises  from  two  facts.  First, 
many  natural  phenomenon  have  been  found  to  be  scaling  over  large  ranges 
of  sizes.  Second,  the  use  of  probability  distributions  to  describe  the 
likelihood  of  observing  specified  differences  in  value  has  proven  to  be 
a  useful  tool  in  the  analysis  of  diffusing  scalars  and  turbulence. 


E.  Relevance  to  Conventional  Fluid  Flow  Modeling 


Typically,  numerical  modeling  of  fluid  flow  involves  processes  that 
will  affect  the  larger  scale  features  of  the  flow,  but  are  too  small  to 
be  resolved  by  the  grid  that  is  used.  Even  the  vortex  methods  (e.g. 
Leonard,  1985)  which  provide  more  detail  in  areas  with  the  most  fine 
structure  face  computational  limitations  that  force  approximations  for 
small  scale  features.  In  general,  the  problem  is  to  define  either  the 
small  scale  features,  or  their  effects,  in  terms  of  those  features  that 
are  resolved  by  the  modeling.  Fractals  have  promise  for  providing 
descriptions  of  small  scale  processes,  but  to  date  there  does  not  seem 
to  have  been  any  practical  application  described  in  the  literature. 
There  have  been  methods  cescribed  for  generating  fields  with  the  appro¬ 
priate  fractal  dimension  (e.g.  Lovejoy  and  Mandelbrot,  1985).  If  fields 
generated  in  this  way  were  realistic  they  could  presumably  provide  a 
basis  for  developing  empirical  parameter izations.  Another  possible 
avenue  can  be  found  in  the  spectral  representations  of  scaling  fields 
(e.g.  Pentland,  1984).  In  principle,  scaling  allows  the  spectral  den¬ 
sity  functions  to  be  extrapolated  to  higher  wave  numbers  (smaller 
scales),  but  it  is  not  clear  whether  computationally  efficient  methods 
can  be  developed  for  taking  advantage  of  this,  or  whether  the  effects  of 
cross-product  terms  (those  involving  more  than  one  parameter)  can  be 
treated  realistically. 


Ill  FRACTALS,  TURBULENCE,  AND  PHYSICAL  SPACE 


A.  Discussion  of  Fractal  Dimension 


1 .  Definitions  and  Hypotheses 

This  section  will  discuss  isotropic  fractals.  Later,  the  concept 
will  be  extended  to  the  anisotropic  case  which  seems  more  applicable  to 
processes  in  the  atmosphere. 

Schertzer  and  Lovejoy  (1983)  show  that  for  many  atmospheric  proper¬ 
ties  (e.g.  kinetic  energy,  the  logarithm  of  the  potential  temperature) , 
the  upper  tail  of  the  probability  distribution  discussed  earlier  falls 


off  algebraically,  rather  than  exponentially.  That  is,  it  can  be 
represented  by: 

Pr( Ac  >  q)  =  k  Ac  q“a  (III-1) 

The  above  relationship  means  that  the  probability  of  a  random  fluctua¬ 
tion  Ac  exceeding  a  fixed  value  q  is  proportional  to  that  fixed  value 
raised  to  the  -a  power,  i.e.  q_a.  The  exponent,  -a  is  a  measure  of  the 
intermittency  of  the  distribution. 

We  can  examine  how  the  moments  behave  as  a  function  of  ot.  Recall¬ 
ing  the  definition  of  the  rth  moment,  vr 

00 

p  =  J  cr  f(c)dc  (II 1—2 ) 

o 

where  f(c)  is  the  probability  density  function,  i.e.  the  derivative  with 
respect  to  c  of  the  cumulative  probability  function  discussed  earlier. 
When  Pr(c  >  q)  =  kq-a,  f(c)  is  given  by 

f (c )  =  kac“(a+1)  ( I I I- 3 ) 

substituting  the  definition  for  the  rth  moment  gives 

°r  -  kaf"  cr“°~1dc  .  (III-4) 

'  o 

A  consequence  of  this  distribution  is  that  there  will  be  a  problem 
in  defining  moments  higher  than  a.  According  to  Schertzer  and  Lovejoy 
(1983),  the  higher  moments  will  diverge  so  that  the  experimental  estima¬ 
tion  will  increase  without  limit  as  the  sample  size  increases. 

There  is  an  interesting  conjecture  that  can  be  made  in  connection 
with  the  divergence  of  the  higher  moments.  I  have  not  seen  the  follow¬ 
ing  line  of  reasoning  spelled  out  in  the  literature,  but  it  is  relevant 
to  the  feasibility  of  the  higher  order  turbulent  closure  schemes  that 
are  included  in  some  fluid  models.  Higher  order  closure  schemes  require 
a  parameterization  of  cross-product  terms  whose  order  may  be  the  same 
as,  or  greater  than  the  order  at  which  the  moments  begin  to  diverge.  It 
seems  reasonable  to  assume  that  if  the  "self-products"  involved  in 
calculating  the  moments  do  not  converge,  there  may  be  problems  with  the 
convergence  of  the  cross-products.  As  noted,  this  is  purely  a  conjec¬ 
ture,  but  it  seems  prudent  to  investigate  it  further  before  elaborate 
higher-order  turbulent  closures  are  developed. 
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2.  Methods  for  Determining  Fractal  Dimension 

Several  different  approaches  have  been  developed  for  calculating 
fractal  dimension.  One  of  these  follows  more  or  less  directly  from  the 
fact  that  fractal  dimension  is  a  measure  of  the  space-filling  properties 
of  a  curve  or  a  higher  dimensional  surface.  Another  approach  estimates 
the  parameters  of  the  probability  distributions,  which  are  directly 
related  to  the  fractal  dimension.  The  fractal  dimension  can  also  be 
derived  from  the  spectral  density  function  and  one  method  is  based  on 
this  approach.  Finally,  it  is  possible  to  estimate  the  fractal  dimen¬ 
sion  of  the  attractor  for  a  given  system  from  a  time  series  at  a  single 
point.  The  applicability  of  this  latter  method  to  fluid  flow  modeling 
applications  is  not  clear,  but  for  the  sake  of  completeness  a  brief 
description  of  the  approach  is  included. 


a.  Subcell  counting  methods 


Mandelbrot  (1975)  used  a  discussion  like  that  which  follows  to 
define  fractal  dimension  and  its  relation  to  the  "space  filling"  con¬ 
cept.  If  the  distribution  of  some  parameter  is  continuous,  then  there 
are  iso-surfaces  containing  all  the  points  that  have  the  same  value.  In 
order  to  determine  the  fractal  dimension  of  one  of  these  iso-surfaces, 
we  select  a  large  cube  in  space  that  contains  at  least  a  part  of  that 
surface.  We  then  define  the  length  of  the  side  of  that  cube  as  an 
external  scale,  L,  and  subdivide  the  cube  into  smaller  cells  with  sides 
of  length,  l.  We  then  count  the  number  of  these  smaller  cells  that 
contain  at  least  one  point  on  the  iso-surface.  If  our  "surface"  is  a 
straight  line,  then  the  number  of  small  cells  through  which  it  passes 
will  be  proportional  to  (L/A) 1 .  If  the  surface  is  a  flat  plane,  then 
the  number  of  small  cells  through  which  it  passes  will  be  approximately 
proportional  to  (L/&)  .  Finally,  if  the  distribution  of  the  parameter 
in  space  is  such  that  every  point  has  the  same  value,  then  the  "iso¬ 
surface"  will  be  solid  and  the  number  of  cells  intercepted  will  be 
proportional  to  (L/ftP:  In  all  these  examples,  the  exponent  of  the  term 
( L/ A)  is  simply  the  Euclidian  dimension  of  the  "surface." 


t 
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The  concept  of  dimension  outlined  in  the  preceding  paragraph  can  be 
extended  to  fractional  dimensions,  if  we  imagine  a  line  with  many 
"wiggles",  and  with  wiggles  upon  those  wiggles  and  so  on  and  on.  In 
such  a  case,  the  number  of  small  cells  (n)  through  which  the  line  passes 
will  be  approximated  by 

n  =  k  (L/A)D  ,  (III-5) 


where  D  is  a  number  larger  than  one,  but  less  than  two.  Similarly,  a 
plane  might  have  roughness  elements  with  smaller  roughness  elements  upon 
them  and  so  on,  in  that  case,  the  exponent  D  would  be  larger  than  two 
but  less  than  three.  The  exponent  D  in  these  cases  is  analogous  to  the 
Euclidian  dimension,  except  that  it  need  not  be  an  integer.  Mandelbrot 
coined  the  term  fractal  (from  "fractional  dimension")  for  geometric 
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shapes  (lines,  surfaces,  etc.)  that  have  this  property.  Fractals  are 
self-similar  in  that  the  fluctuations  at  a  small  scale  are  related  to 
those  at  a  larger  scale  by  a  constant  power  of  the  scale  ratio  over  a 
wide  range  of  scales.  Love joy  and  Schertzer  (1986)  have  pointed  out 
that  there  are  good  physical  reasons  for  the  exponent  for  vertically 
oriented  features  to  differ  from  that  from  those  for  the  horizontal 
features  in  the  atmosphere.  This  will  be  discussed  further  in  a  subse¬ 
quent  section. 

According  to  Greenside  et  al  (1982),  Takens  (1981),  suggested  an 
algorithm  for  computing  fractal  dimension  based  on  Mandelbrot's  (1975) 
definition  outlined  above.  Greenside  et  al  (1982)  concluded  that  this 
approach  was  impractical  for  systems  with  dimensions  greater  than  2, 
because  of  convergence  problems.  These  findings  do  not  seem  to  preclude 
the  use  of  the  method  for  finding  the  fractal  dimension  for  isopleths  of 
a  scalar  in  a  plane,  although  the  convergence  may  be  slow.  The  method 
should  probably  not  be  applied  to  problems  involving  high-dimensional 
phase  spaces.  When  the  definition  is  reduced  to  two  dimensions,  a 
square  and  smaller  squares  replace  the  cube  and  subcubes,  so  it  can  be 
used  to  examine  isopleths  of  a  scalar  in  a  plane  intersecting  an  iso¬ 
surface. 

The  cell  counting  concept  is  reasonably  direct,  and  easily  applied 
to  an  array  of  discrete  numbers  such  as  a  digitized  image.  This  makes 
it  particularly  useful  for  the  analysis  of  cloud  photos,  or  backscatter 
cross-sections  through  smoke  plumes,  as  are  obtained  with  mobile  lidar 
systems.  Another  potential  application  would  be  to  digitized  images  of 
the  markers  in  laboratory  turbulence  experiments.  One  interesting 
feature  of  the  technique,  as  it  is  applied  to  smoke  plume  cross-sections 
arises  from  the  fact  that  the  isopleths  tend  to  be  concentric  (Ludwig 
and  Nitz,  1986),  so  the  fractal  dimension  obtained  for  low-value  iso¬ 
pleths  are  applicable  at  the  edges  of  the  plume,  while  those  for  the 
higher  values  will  be  characteristic  of  the  center. 

Personal  experience  has  shown  that  one  does  not  have  complete 
freedom  in  the  choice  of  the  area  over  which  an  algorithm  based  on  the 
cell  counting  method  can  be  applied,  especially  in  the  case  of  arrays  of 
discrete  values.  The  method  uses  linear  regression  to  find  the  con¬ 
stants  for  the  logri;hmic  form  of  Equation  III-5,  i.e. 

log(n)  =  log(k)  +  D  log(L/i.) .  (III-6) 

In  practice,  we  want  as  many  values  of  N  as  possible.  For  a  discrete 
array,  each  N  must  correspond  to  an  integer  value  of  ih/l).  This  means 
that  we  must  select  the  large  squares  so  that  they  have  a  side  length  L 
that  has  many  divisors.  For  example,  if  we  choose  a  large  square  whose 
side  L  is  a  prime  number  of  grid  points,  then  we  will  only  get  two 
points  from  which  to  calculate  the  regression  constants  in  Equation 
1 1 1  -  6 .  It  turns  out  that  a  wide  choice  of  values  can  be  achieved  by 
judicious  use  of  products  of  powers  of  two,  three  and  five.  Using  this 
approach,  it  is  possible  to  limit  the  analysis  to  a  region  which 
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encompasses  that  part  of  the  image  which  is  of  interest,  but  does  not 
extend  very  far  into  background  areas. 


b.  A  Method  Based  on  Spatial  Variation  Probability 
Distributions 

We  have  already  noted  that  spatially  varying  fields  are  frequently 
characterized  by  the  probability  distributions  for  the  difference 
between  values  at  two  points  separated  by  a  specified  distance,  and  that 
there  is  a  functional  relationship  for  the  probability  distributions  for 
different  separation  distances.  Schertzer  and  Love joy  (1983)  presented 
the  following  approach  that  makes  use  of  these  definitions  to  determine 
fractal  dimension.  Starting  with  the  definition  of  a  fractal  (scaling) 
pattern,  i.e. 

Pr[Ac (AY)  >  q]  =  Pr[AHAc(Y)  >  q]  (III-7) 

As  before,  Ac(x)  is  a  change  in  value  of  the  variable  c  over  the  dis¬ 
tance  x.  In  practice,  Y  will  generally  be  the  smallest  separation 
between  points  in  the  available  data  set.  Equation  (III-7)  can  be  used 
to  estimate  the  exponent  H  if  one  determines  various  quantile  values  as 
a  function  of  the  separation  distance  (AY)  and  then  determines  the  best 
fit  value  for  H. 

If  we  accept  the  hypothesis  of  Schertzer  and  Love joy  (1983)  that 
the  upper  tail  of  the  distribution  has  a  hyperbolic  form  for  large  q 

Pr (Ac  >  q)  =  kq"a  ,  (III-8) 

then  the  constants  H  and  a  can  be  estimated  from  the  data  as  the  best 
fit  values  for  the  following  relationship — again,  for  large  values  of  q: 

Pr[ Ac ( AY )  >  q]  -  k  AHq“a  (III-9) 

The  steps  necessary  to  apply  this  method  are  as  follows: 

(1)  Determine  AH  and  a  at  different  separations  for  upper  tail 
percentiles  less  than  a  few  percent. 

(2)  Determine  the  average  value  of  a,  providing  that  the  indivi¬ 
dual  values  of  a  are  sufficiently  similar. 

U 

(3)  Find  H  by  a  log-linear  regression  of  An  values. 

The  scaling  parameter  H  is  related  to  fractal  dimension  as  follows 
(Pentland,  1984); 

D  ■=  2  -  H,  ( I II- 1 0 ) 

Schertzer  and  Lovejoy  (1983)  analyzed  distributions  of  several  different 
atmospheric  parameters  by  the  above  method  and  found  that  the  data 
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support  the  hypothesis  of  a  hyperbolic  upper  tail  for  the  probability 
distribution.  They  also  found  that  the  data  showed  that  the  distri¬ 
butions  are  anisotropic  with  vertical  fractal  dimensions  that  are  dif¬ 
ferent  from  those  in  the  horizontal.  The  importance  of  this  fact  will 
be  discussed  later. 


c .  Spectral  Analysis  Methods 


The  method  described  in  the  preceding  section,  based  on  the  proba¬ 
bilities  of  finding  differences  of  a  specified  magnitude  as  a  function 
of  separation  distance,  is  closely  related  to  the  estimation  of  spatial 
autocorrelation.  Therefore,  we  might  expect  that  there  would  be  spec¬ 
tral  analysis  methods  that  could  be  used  to  estimate  fractal  dimensions. 
For  a  variable  that  satisfies  the  scaling  relationship  presented  above, 
the  exponent  of  the  corresponding  power  spectrum  (-6)  is  related  to  the 
scaling  exponent  H  as  follows: 

6  =  2H  +  1  .  (III-11) 
Combining  Equations  III-10  and  II I  —  1 1  gives 

B  *  5  -  2D  .  (III-12) 

Thus,  the  power  spectrum  P(f)  of  a  scaling  field  will  be  proportional  to 
f2D  ,  where  P  is  the  spectral  density  as  a  function  of  wave  number  f. 
This  relationship  can  be  used  as  the  basis  of  a  method  for  estimating 
the  fractal  dimension  of  a  field,  where  fast  Fourier  Transform  (FFT) 
methods  are  used  to  generate  the  amplitudes  of  the  Fourier  terms  in  wave 
number  space.  The  power  spectrum  is  then  fit  without  regard  to  direc¬ 
tion.  In  essence,  the  procedure  averages  the  amplitudes  around  circles 
of  constant  wave  number  radius.  If  anisotropic  effects  are  to  be 
included,  the  power  spectra  along  different  axes  would  be  examined. 
Pentland  (198*1)  has  shown  that  the  spectral  method  can  be  applied  to 
fairly  small  subsections  of  an  image.  He  used  the  method  to  charac¬ 
terize  the  fractal  dimensions  of  16  x  16  pixel  subsections  of  photo¬ 
graphs.  The  fractal  analysis  provided  a  basis  for  distinguishing  the 
inhomogeneities  in  the  areas  of  the  picture  that  had  similar  texture. 
The  method  is  efficient  because  it  uses  FFT  algorithms  and  simple  linear 
regression,  but  substantial  changes  would  be  needed  to  address  aniso¬ 
tropic  effects. 


d.  Estimation  of  Fractal  Dimensions  from  Time  Series 


As  we  noted  earlier,  the  time  evolution  of  some  system  can  be  fully 
described  by  a  trajectory  in  an  n-dimensional  phase  space.  Trajectories 
for  atmospheric  systems  are  often  contained  within  "surfaces"  whose 
dimension  d  is  smaller  than  that  of  the  complete  phase  space,  i.e. 
d  <  n.  Measurements  are  more  commonly  available  for  a  long  period  of 
time  at  a  few  individual  locations  than  at  enough  points  to  describe  the 
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state  of  the  system.  Therefore,  it  is  be  desirable  to  have  a  method 
that  can  use  time  series  data  to  estimate  fractal  dimension  of  the 
surface  containing  the  state  trajectories.  Fraedrich  (1986)  presents  a 
discussion  of  how  this  is  done.  This  discussion  and  that  of  Packard  et 
al  (1980)  form  the  basis  for  the  following  description.  If  the  dynamics 
of  the  system  are  described  by  a  set  of  n  ordinary  differential  equa¬ 
tions 


XJ  ’  W  *2 . J  ■  ' . " 


(III- 1 3) 


then  the  phase  space  is  described  by  the  n  coordinates  Xj  ,  (j  =  1,  ..., 
n).  These  are  the  dependent  dynamical  variables.  The  symbol  (*) 


indicates  differentiation  with  respect  to  time;  the  time  evolution  of 
the  system  can  be  written  as  an  n-dimensional  vector  function  of  time 
[x.|(t) . xn(t)j,  defining  the  trajectory  in  phase  space.  The  dif¬ 
ferential  Equations  (III— 1 3)  can  be  reduced  to  a  single  (albeit  highly 
nonlinear)  differential  equation  involving  only  one  of  the  variables,  if 
all  the  others  are  eliminated  by  differentiation.  This  nth  order  equa¬ 
tion  will  be  of  the  form 


:(n)  =  f[x,  x' . x^1)] 


(111-14) 


This  is  equivalent  to  n  equations  describing  the  time  series  of  the 
variable  x(t)  and  its  first  n-1  derivatives,  x'(t),  ...,  x^n_1^(t).  The 
Equations  (III-l 4)  are  equivalent  to  the  original  set  of  differential 
Equations  (III-13),  hut  with  the  important  difference  that  all  the 
variables  can  be  determined  from  measurements  at  a  single  point  —  if 
the  measurements  have  sufficient  temporal  resolution  to  obtain  the 
necessary  higher-order  derivatives. 


According  to  Fraedrich  (1986),  if  the  dimension  d  of  the  attractor 
contained  within  the  original  n-dimensional  phase  space  is  much  less 
than  n,  it  can  then  be  described  in  the  new  phase  space  made  up  from  the 
single  variable  and  its  derivatives.  In  fact,  Fraedrich  (1986)  notes 
that  the  dimensionality  of  the  new  phase  space  can  be  smaller  than  that 
of  the  original  phase  space,  i.e.  one  need  not  include  all  n-1  deriva¬ 
tives.  In  fact,  2d  +  1  variables  (i.e.  the  variable  and  its  first  2d 
derivatives  should  be  sufficient). 


In  practice,  one  does  not  need  to  know  the  number  of  state  vari¬ 
ables  in  the  original  phase  space,  if  sufficiently  many  derivatives  are 
included.  Also,  in  practice  the  measurements  are  available  at  discrete 
times,  so  that  (2d-1)  discrete  time  lag  multiples  are  substituted  for 
continuous  derivatives.  In  this  new  system,  we  deal  with  the  following 
vector 


2(t)  =  (x(t),  x(t+i),  ...,  x [t  +  ( 2d- 1 ) t ] } 


( III- 1 5 ) 


In  order  to  ensure  linear  independence,  the  time  shift  t  is  chosen  so 
that  autocorrelations  vanish. 
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Observed  time  series  of  a  single  variable  describe  trajectories  in 
this  new  phase  space,  so  that  one  of  the  methods  discussed  earlier  can 
be  used  to  determine  the  dimension  of  the  structure  formed  by  those 
trajectories.  Of  course,  one  must  know  that  dimension  first  to  deter¬ 
mine  how  many  derivatives  (time  shifts)  must  be  included  in  the  calcula¬ 
tions.  It  is  desirable  to  minimize  that  number  because  the  required 
calculations  increased  dramatically  with  increasing  numbers  of  dimen¬ 
sions.  Also,  as  noted  before,  the  cell  counting  method  may  not  be 
appropriate  for  dimensions  greater  than  two  (Greenside  et  al,  1982). 

Fraedrich  (1986)  uses  an  algorithm  of  Grassberger  and  Procaccia 
(1983)  that  provides  an  estimate  for  the  lower  limit  of  the  fractal 
dimension.  Fraedrich  (1986)  notes  that  in  most  of  the  cases  that  were 
analyzed,  this  lower  limit  was  very  close  to  the  actual  value.  In  order 
to  circumvent  the  problems  arising  from  a  lack  of  foreknowledge  about 
the  fractal  dimension,  Fraedrich  (1986)  calculates  fractal  dimension 
using  successively  larger  numbers  of  derivatives  until  he  finds  that  the 
calculated  dimension  no  longer  changes  with  an  increase  in  the  number  of 
dimensions  in  the  phase  space. 

Obviously,  methods  based  on  time  series  from  a  single-point  have 
considerable  practical  advantage  over  those  that  must  analyze  data  from 
a  multitude  of  individual  points.  Thus,  this  method  deserves  further 
exploration  with  regard  to  its  application  to  studies  of  turbulence, 
including  results  obtained  from  numerical  simulations,  especially  large 
eddy  simulations  (LES).  To  date,  it  appears  that  the  applications  of 
the  technique  have  been  much  like  those  of  Fraedrich  (1986)  and  have 
involved  climatic  and  other  variabilities  on  time  scales  much  longer 
than  those  appropriate  to  studies  of  turbulence  and  turbulent  effects. 


3.  Generation  of  Fractal  Fields  with  Specified  Characteristics 

One  possible,  practical  use  of  the  fractal  concept  would  be  to 
develop  parameter izations  for  estimating  fractal  characteristics  from 
large-scale  flow  parameters  so  that  fields  with  the  appropriate  charac¬ 
teristics  could  be  generated  and  used  for  empirical  characterization  of 
terms  appearing  in  conventional  closure  schemes.  This  approach  would 
require  a  method  for  generating  fields  with  the  appropriate  fractal 
characteristics.  This  problem  has  been  addressed,  in  large  part  because 
of  its  relevance  to  the  generation  of  realistic  looking  artificial 
landscapes . 

Probably  the  most  efficient  way  of  generating  such  fields  is  to  use 
the  appropriate  Fourier  series,  but  that  technique  may  not  allow  intro¬ 
duction  of  randomness.  Lovejoy  and  Mandelbrot  (1985)  describe  a  tech¬ 
nique  that  they  call  the  "fractal  sums  of  pulses"  (FSP)  method  for 
generating  rainfall  intensity  statistics  that  have  the  appropriate 
fractal  properties.  In  one  dimension,  i.e.  a  time  series  of  rainfall 
intensity,  R(t),  the  FSP  method  works  as  follows.  The  function  R(t)  is 
generated  as  a  sum  of  randomly  located  rectangular  pulses  with  random 


height  (rainfall  intensity)  and  duration.  The  random  locations  of  the 
centers  are  distributed  according  to  a  Poisson  process  with  constant 
rate,  v.  To  ensure  scaling,  Lovejoy  and  Mandelbrot  (1985)  choose 
Pr(p'  >  p)  to  be  proportional  to  p'  ;  p  is  the  duration.  Thus,  the 
probability  of  a  random  duration  p'  exceeding  p  is  proportional  p-1. 
The  height  of  the  rectangle  (or  rainfall  intensity  increment,  AR)  is 
taken  to  be  equal  to  ±  (p  'y^a.  They  point  out  that  for  any  given 
duration  p' ,  the  pulse  intensity  has  a  random  sign,  but  a  fixed  absolute 
value. 

This  one  dimensional  model  scales  with  an  exponent  H=1/a.  A  change 
in  duration  scale  by  a  factor  of  X  will  change  the  intensity  by  a  factor 
of  A 1  ^a.  The  intensity  probability  distribution  will  be  hyperbolic  with 
an  exponent  of  a,  because  of  the  rainfall  duration  probability  distribu¬ 
tion  and  choice  of  relationships  between  rainfall  duration  and  rainfall 
intensity . 

The  pulses  need  not  be  rectangular,  but  can  be  smooth  and  continu¬ 
ous  if  their  length  scales  as  A  and  the  intensity  scales  as  A  /ct. 
Lovejoy  and  Mandelbrot  (1985)  suggest  a  pulse  shape  with  an  amplitude  AR 
proportional  to  AR  exp[- (u/p) ,  where  u  is  the  distance  from  the  pulse 
center.  This  arbitrary  shape  is  convenient  because  the  pulse  is  smooth 
for  small  values  of  S,  becoming  more  discontinuous  (more  nearly 
rectangular  pulses)  as  S  becomes  larger. 

The  method  can  be  extended  to  more  dimensions.  For  an  isotropic 
distribution  in  two  dimensions,  the  simplest  pulse  is  a  randomly  placed 
right  circular  cylinder.  The  probability  that  the  base  area  A  of  the 
cylinder  exceeds  a  specified  value  is  inversely  proportional  to  the 
value,  while  cylinder  height  (intensity)  equals  ±  A1  .  It  is  necessary 
to  define  an  intensity  threshold  and  set  values  below  that  threshold  to 
zero  when  a  positive  field  (e.g.  rainfall)  is  represented. 

Although  the  rain  fields  constructed  by  the  above  method  have  the 
generally  appropriate  fractal  properties,  they  tend  to  be  sharp-edged 
and  produce  too  few  separate  rain  areas.  Lovejoy  and  Mandelbrot  (1985) 
go  on  to  describe  other  models  using  pulses  with  annular  bases  and 
elliptical  annuli  of  varying  eccentricity  to  simulate  anisotropic 
effects. 

Medler  et  al  (1986)  have  used  some  of  the  ideas  presented  by 
Lovejoy  and  Mandelbrot  to  generate  reasonably  realistic  images  of  smoke 
plumes.  However,  it  should  be  noted  that  their  simulations  are  not 
based  on  observed  fractal  properties.  They  have  simply  chosen  fractal 
properties  to  achieve  a  degree  of  visual  realism  in  the  images  that  they 
produce.  Lovejoy  and  his  coworkers  (e.g.  Lovejoy  and  Schertzer,  1985, 
1986;  Lovejoy  and  Mandelbrot,  1985;  Schertzer  and  Lovejoy,  1983)  have 
produced  images  that  simulate  clouds  and  precipitation  areas  using  the 
methods  described  above.  Unlike  Medler  et  al  (1986),  they  have 
attempted  to  use  parameters  based  on  observational  evidence.  They  have 
also  attempted  to  include  anisotropic  effects. 
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Although  considerable  progress  has  been  made  in  producing  scalar 
fields  that  have  fractal  properties  similar  to  those  observed  for 
certain  atmospheric  variables,  there  are  some  very  important  problems 
which  do  not  yet  appear  to  have  been  addressed.  In  particular,  I  have 
found  no  reports  where  fractals  have  been  used  extend  to  some  larger 
scale  description  to  smaller  scales.  Medler  et  al  (1986)  and  Pentland 
(personal  communication,  1986)  have  superimposed  fractal  patterns  on 
general  Gaussian  plume  concentration  cross  sections,  but  have  not  used 
observed  fractal  parameters  in  the  process.  Furthermore,  both 
applications  have  been  limited  to  isotropic  fractals. 


Another  area  where  much  remains  to  be  done  with  fractals  is  vector 
field  simulation.  This  represents  a  formidable  task,  but  is  very  impor¬ 
tant  because  appropriate  vector  fields  could  serve  as  a  basis  for  turbu¬ 
lent  closure  or  direct  simulation  of  subscale  eddy  processes.  Although 
at  least  one  method  for  estimating  fractal  dimension  can  be  applied  to 
vector  quantities  (e.g.  Pentland,  1984),  it  is  not  immediately  obvious 
how  that  definition  can  be  used  to  generate  vector  fields  with  appropri¬ 
ate  fractal  properties  and  statistical  relationships  between  components. 


B.  Why  Turbulence  Might Be  Expected  to  be  Well  Described  by  Fractals 


To  this  point,  the  discussion  has  focused  on  definitions  of  frac¬ 
tals  and  related  concepts  without  giving  serious  attention  to  the  physi¬ 
cal  reasons  why  they  may  be  important  in  the  study  of  turbulence  and 
turbulent  diffusion  of  scalars.  There  are  at  least  five  kinds  of  evi¬ 
dence  which  suggest  the  connection: 


(1  ) 


(2) 

(3) 

(4) 

(5) 


The  structure  of  turbulent  eddies  observed  in  flow  visualiza¬ 
tion  experiments 

Similarity  arguments  from  classical  turbulence  theory 
Detailed  numerical  flow  simulations 
Characteristics  of  the  governing  equations,  and 


detailed  analysis  of  the  spatial  distribution  of  atmospheric 
scalars. 


These  are  discussed  in  the  following  sections, 


1 


Shear  Layer  Flow  Visualization 


Reynolds  (1985)  provides  a  qualitative  discussion  of  the  processes 
and  the  structure  in  a  turbulent  mixing  layer  forming  in  the  shear  zone 
between  fluids  flowing  at  different  velocities.  The  structure  evident 
in  Figure  1  occurs  quite  commonly.  Although  Reynolds  does  not  mention 
fractals,  his  discussion  provides  evidence  of  scaling  and  similarities 
between  the  structure  of  the  motion  field  on  different  scales.  What 
follows  is  based  on  the  material  presented  by  Reynolds  (1985). 


Figure  D-1  shows  what  happens  in  a  mixing  layer  flow.  Initially, 
vorticity  flows  from  the  tip  of  the  separator  that  is  used  to  generate 
the  flow.  An  unstable  2-dimensional  shear  layer  is  formed.  The  insta¬ 
bility  can  be  excited  by  random  noise  or  slight  vibrations.  Once 
excited,  it  grows  rapidly  so  that  vorticity  is  quickly  concentrated  in 
nearly  discrete  vortices,  oriented  generally  perpendicular  to  the  flow 
direction,  although  they  may  include  some  irregularities .  It  is  impoi — 
tant  to  note  that  not  all  the  vorticity  is  concentrated  in  the  major 
features,  but  some  remains  in  "braids"  that  connect  the  major  vortices, 
as  shown  in  Figure  D-1.  As  will  become  evident,  this  has  important 
consequences . 

As  noted  above,  there  may  be  irregularities  in  the  spacing  or 
orientation  of  the  main  vortices.  Such  irregularities  can  cause  the 
vortices  to  move  at  different  velocities,  so  that  they  merge  and  their 
size  increases  in  the  downstream  direction.  Figure  D-2  (adapted  from 
Reynolds,  1985)  shows  schematically  how  this  takes  place.  The  regularly 
spaced,  uniform  vortices  at  the  left  of  the  figure  interact  so  that 
induced  motions  (e.g.  the  motion  of  vortex  C  caused  by  vortex  B)  are 
cancelled  (e.g.  by  vortex  D).  Some  of  the  regularly  spaced  vortices 
shown  at  the  right  of  the  figure  will  move  up  and  others  will  move 
down.  Eventually  some  pairs  move  close  together  and  orbit  one  another 
or  merge . 

It  is  in  the  braids  between  the  large  scale  vortices  where  the 
evidence  of  fractal  structure  may  be  found.  As  noted  before,  some 
vorticity  remains  in  the  braids  to  be  stretched  as  the  large  scale 
vortices  "wind  in"  the  fluid  between  them.  Stretching  intensifies  the 
vorticity  in  the  braids,  forming  new  vortices  with  axes  aligned  along 
the  principal  strain  direction  as  shown  in  Figure  D-3.  These  vortices 
undergo  the  same  processes  described  above,  but  on  a  smaller  scale  and 
with  different  orientation.  One  can  expect  that  the  new,  smaller  vor¬ 
tices  also  contain  irregularities  that  cause  them  to  merge,  and  that  new 
"minibraids"  form  and  are  stretched  along  new  strain  axes.  These  new 
structures  can  undergo  similar  deformation  on  a  still  smaller  scale,  so 
that  the  processes  are  likely  to  continue  down  to  the  scale  of  viscous 
dissipation.  Thus,  qualitatively  we  would  expect  scaling  and  a  similar¬ 
ity  of  turbulent  structure  to  extend  over  a  wide  range  of  scales,  from 
the  outer  scale  defined  by  the  large,  merged  vortices  down  to  the  scale 
of  molecular  dissipation  processes.  It  should  be  noted  that  analogous 
reasoning  can  be  applied  to  other  flow  types  such  as  jets,  boundary 
layers  and  wakes  which  are  also  characterized  by  regions  of  strong  shear 
and  concentrated  vorticity. 


2.  Classical  Theory 

The  term  "classical  theory"  in  this  section's  title  refers  to 
Kolmogorov's  early  deductions  concerning  the  cascade  of  energy  through 
the  turbulent  spectrum.  However,  the  derivation  that  follows  is  taken 
from  Frisch  et  al  (1978),  because  their  discussion  provides  a  somewhat 
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more  physical  picture  of  the  process.  Similar  arguments  have  also  been 
given  by  Love joy  and  Schertzer  (1986).  We  begin  by  defining  the  energy 
spectrum  E(k)  as  the  kinetic  energy  per  unit  mass  per  unit  wave  number 
k.  For  purposes  of  argument,  Frisch  et  al  (1978)  choose  to  define  a 
discrete  spectrum  of  eddies  beginning  with  the  large  scale  £,Q  where  the 
energy  is  introduced  and  proceeding  to  successively  smaller  wave  numbers 
according  to  the  relationship 

i  =  l  2~n  ,  n-0,1,2,  ...,  (III-16) 

n  o 

the  corresponding  wave  numbers  are 

kn=1/i,n  .  (Ill— 17) 

The  discretized  kinetic  energy  per  unit  mass  in  the  wave  numbers  near  Jln 
is  given  by 


E 


n 


E(k )dk 


(III— 18) 


If  the  turbulence  is  statistically  stationary,  energy  introduced  at 
large  scales  ( 2,Q )  is  transferred  to  successively  smaller  scales  until 
the  dissipation  scale  J.d  is  reached.  We  can  define  a  characteristic 
velocity  vn  in  terms  of  En  for  the  kn  eddies,  i.e. 

E  -  v  2  .  (Ill— 19) 
n  n 

It  should  be  noted  that  vn  is  not  total  velocity,  but  rather  a  velocity 
difference  characteristic  of  that  particular  eddy  size  i,n.  With  this  in 
mind,  an  eddy  turnover  time  tn  can  be  defined 


t  -  l/vn  .  (III-20) 
n  n  n 

Frisch  et  al  (1978)  go  on  to  argue  that  in  the  inertial  subrange — where 
tn  is  much  greater  than  the  viscous  dissipation  time,  &n  'V  (v  is  the 
kinematic  viscosity),  and  much  less  than  the  characteristic  time  for  the 
larger  scale  motions,  £0/v0 — we  can  define  an  energy  (per  unit  mass) 
transfer  rate  from  eddies  of  wave  number  kn  to  kn+i  to  be 


e  -  E  /t  -  v  /% 
n  n  n  n  n 


(II 1—2 1 ) 


If  the  process  is  stationary  as_hypothesized ,  then  the  energy  dissipa¬ 
tion  rate  en  will  be  a  constant  e,  and  we  can  now  solve  for  vn  and  En: 


-  (a  )1/3 

(III-22 ) 

n 

-  «*//3 

9 

(III-23) 
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Equations  (III-22)  and  C I I I— 23 )  are  the  same  as  Kolmogorov's  result  for 
the  structure  functions.  A  Fourier  transformation  provides  the  wave 
number  spectrum 


E(k )  -  e  2/3  k“5/3 


We  can  also  obtain  the  eddy  turnover  time 


t  -  e 
n 


(111-24) 


(HI-25) 


The  above  results  can  be  rewritten  to  show  their  scaling  properties 
more  clearly.  For  example  from  Equation  (III-22) 


E  -  e  2/3  £  2/3  -  e  2/3  £  2/3  (£  /£  )2/3 
n  n  o  n  o 


E  -  E[(£  /£  )  £  ]  -  E  ( £/£)“ 
n  n  o  o  o  n  o 


(III-26) 


(III-27) 


Thus,  the  energy  per  unit  mass  per  unit  volume  scales  according  to 
(£/£q)2/3.  There  is  a  tacit,  but  very  important  assumption  in  the  above 
relationship,  i.e.  all  eddy  sizes  are  assumed  to  be  spread  more-or-less 
uniformly  throughout  the  same  volume. 

For  greater  generality,  Frisch  et  al .  (1978)  assume  that  the 
smaller  eddies  are  less  space  filling  than  the  larger  ones.  Certainly, 
the  qualitative  behavior  of  mixing  layer  turbulence  described  in  the 
preceding  section  makes  such  an  assumption  plausible.  This  assumption 
requires  that  we  introduce  a  new  parameter  to  characterize  the  degree  to 
which  the  cascade  is  space  filling.  If  the  whole  volume  is  filled  by 
eddies  of  all  sizes,  then  the  cascade  we  have  been  discussing  (where 
£n  «  £~n)  has  2^  times  as  many  eddies  of  size  £n+-|  as  it  has  of  size 
£n.  That  is,  it  takes  eight  eddies  of  dimension  £/2  to  fill  the  volume 
occupied  by  one  of  dimension  £.  Frisch  et  al  (1978)  define  8  to  be  the 
ratio  of  the  average  number  of  £n+1  eddies  produced  by  each  £n  eddy  and 
the  maximum  possible,  i.e.  eight.  If  eddies  of  size  £Q  are  space  fill¬ 
ing,  then  eddies  of  size  £n  will  fill  only  a  fraction  Bn  of  the  total 
space 


8n  «  Bn  =  (N/23)n 


(III-28) 


where  N  is  the  average  number  of  eddies  formed  by  each  eddy  of  the 
preceding  generation;  N  <,  8. 

Frisch  et  al  (1978)  assume  that  the  eddies  of  the  (n+1)th  genera¬ 
tion  are  positionally  correlated  with  those  of  the  nth  generation  by 
embedding  or  attachment.  This  seems  to  mean  that  they  believe  that  the 
smaller  eddies  are  produced  in  the  same  general  area  as  the  larger  eddy 
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from  which  they  arise.  The  discussion  of  mixing  layer  turbulence  given 
in  the  preceeding  section  suggests  that  the  "positional  correlation"  is 
more  likely  to  be  one  of  "attachment"  than  "embedding."  Frisch  et  al 
(1978)  are  assuming  that  the  region  where  an  eddy  is  formed  becomes  an 
"active  region"  for  the  cascade  to  smaller  sizes.  There  seems  to  be  one 
other  tacit  assumption  in  the  arguments  of  Frisch  et  al  (1978)  that  they 
do  not  explicitly  state.  It  appears  that  they  assume  that  the  average 
number  of  eddies  formed,  N,  does  not  vary  systematically  according  to 
eddy  size. 

Frisch  et  al  (1978)  redefine  the  typical  velocity  difference  vn  in 
terms  of  active  regions  only,  so  that  the  kinetic  energy  per  unit  mass 
associated  with  scales  on  the  order  of  £n  is  given  by 


0  v 
n  n 


(II 1—29 ) 


The  characteristic  time  for  energy  transfer  from  £.n~scale  eddies  to 
smaller  scales  remains  the  same,  i.e.  tn  »  in/vn,  assuming  that  the 


assuming  that  the 


production  of  eddies  of  the  next  smaller  scale  arises  from  the  internal 
dynamics  of  the  larger  eddies  producing  them.  This  seems  reasonable  in 
light  of  the  mixing  layer  discussion  given  earlier,  because  smaller 
scale  eddy  generation  depends  on  the  stretching  of  the  braids  caused  by 
the  "winding  up"  in  the  larger  eddies. 

If,  as  before,  we  assume  a  steady-state  condition  so  that  the  rate 
of  energy  transfer  is  independent  of  scale  over  the  inertial  range,  i.e. 


e  -E/t  -  0  v  /£  -  e  , 

n  n  n  n  n  n 


then,  the  following  relationships  are  obtained 


v  -  c-,/3  l1'3  U  /»  )-(3-D)/3 
n  n  n  o 


t  -  r1/3 12/3  u/t)(3-D)/3 

n  n  n  o 


E  -  £2/3  i2/3  d  n  )<3-D,/3 
n  n  n  o 


E<k>  -  :2/3  k‘5/3  <ki  >-<3-°>'3 
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(III-30) 


(III-31 ) 


(III— 32) 


(III-33) 
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In  the  above  expressions,  D  is  defined  by  N=2  .  The  scaling  nature  of 
the  relationships  becomes  evident  if  we  rewrite  Equation  (III —33 ) 

,5  D  5  D 

( - _)  ( _ X 

-2/3  2/3  33  33^  (III— 35) 

E[ ( l  /£)*]-  e  *  l  5  U  /l  )  5  5  ~  U  /l  )  5  5  E 

noo  ono  no  o 
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This  is  the  Kolmogorov  result,  but  with  a  correction  of  D/3  in  the 
exponent  to  account  for  intermittency  in  the  energy  transfer  to  smaller 
scales.  Finally,  the  relationship  among  N,B  and  D  is: 

Bn  -  (N/23)n  =  (2D~3)n  .  (III-36) 

The  above  derivation  demonstrates  the  scaling  properties  of  impor¬ 
tant  turbulent  parameters.  It  also  relates  them  to  classical  theory  and 
some  of  the  observed  physical  factors  that  lead  to  intermittency.  In 
the  process,  the  fractal  dimension  is  introduced  to  describe  the  space¬ 
filling  properties  of  the  "surfaces'’  within  which  transfers  of  energy 
from  larger  to  smaller  scales  take  place. 

It  should  be  noted  that  there  is  a  scale  at  which  the  viscous 

dissipation  takes  place,  and  it  is  different  from  that  given  by 

Kolmogorov.  The  outer  scale  where  energy  is  introduced  to  the  cascade 

2 

remains  the  same,  and  defines  e  -  v  /£  .  At  the  other  end,  we  can 

o  o 

still  define  the  scale  at  which  the  viscous  dissipation  takes  place  by 

equating  turnover  time,  t  -  e1/3  Jl2/3  (8,  /i  )^3  D^/3  to  the  viscous 

2  n  n  n  o 

diffusion  time,  8T  /v.  Solving  gives 


l 


d 


3/O+D) 


(III-37) 


Thus,  when  the  intermittency  effects  are  included,  the  dissipation  scale 

will  differ  from  the  Kolmogorov  microscale  (v3/e)1/,J4.  They  are  equal 
only  when  D  =  3,  i.e  for  the  space-filling  case. 


3.  Numerical  Simulation  of  Small  Scale  Structure 


Most  fluid  simulations  parameterize  turbulent  effects  and  so  pro¬ 
vide  little  evidence  for  any  fractal  structures  related  to  turbulence. 
However,  Chorin  (1982)  performed  a  very  interesting  numerical  experiment 
that  bears  a  close  relationship  to  the  formation  of  vortices  in  the 
braids  between  eddies  that  was  discussed  earlier.  Chorin  (1982)  con¬ 
sidered  a  straight  vortex  with  a  single  perturbation,  embedded  in  a  3- 
dimensionally  periodic  domain.  The  initial  vortex  ran  from  the  center 
of  the  bottom  of  the  box  to  the  center  of  the  top,  with  a  jog  like  that 
showed  schematically  in  Figure  D-M .  Thus,  he  began  with  a  perturbation 
on  a  structure  very  much  like  the  vortex  filaments  in  the  braids. 

Chorin  (1982)  used  a  vortex  element  simulation  technique,  much  like 
those  described  by  Leonard  (1985).  The  overall  structure  was  simulated 
by  a  number  of  straight  vortex  tube  segments.  After  the  simulation 
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began,  the  vorticity  stretched  very  rapidly  and  increased  the  number  of 
segments  that  were  required  to  approximate  the  distorted  vortex  tube. 
The  rapid  increase  in  the  required  number  of  vortex  elements  limited  the 
time  that  the  calculations  could  be  continued.  Figure  D-5  shows  the 
general  configuration  after  10,  20,  30,  and  40  time  steps,  corresponding 
to  elapsed  times  of  0.65,  0.88,  1.04,  and  1.21,  respectively .  Presum¬ 
ably,  though  not  stated,  the  units  are  seconds. 

Figure  D-5  shows  that  the  general  orientation  of  the  pattern 
remains  vertical  while  becoming  very  complex.  Chorin  (1982)  evaluated 
the  Hausdorff  (fractal)  dimension  of  the  structure  in  Figure  D-5  and 
found  it  to  be  on  the  order  of  2.5.  This  is  consistent  with  values 
suggested  by  Mandelbrot  (1977)  for  the  fractal  dimension  of  turbulent 
structures . 


4.  The  Navler-Stokes  Equations 

The  three  preceding  sections  show  that  there  is  a  wide  variety  of 
evidence  to  suggest  the  general  applicability  of  fractal  concepts  to  the 
description  of  turbulence.  Before  proceeding  to  discuss  some  evidence 
gathered  from  atmospheric  observations,  it  is  worth  noting  that  solution 
spaces  for  the  Navier-Stokes  equations  are  widely  believed  to  be 
fractals  as  well.  The  focus  of  this  review  is  not  on  that  aspect  of  the 
problem,  so  it  will  not  be  reviewed  in  detail. 

Lorenz  (1963)  devised  a  set  of  very  simple  differential  equations 
whose  trajectories  in  state  space  are  contained  in  convoluted  surfaces 
that  seem  to  have  many  of  the  same  characteristics  that  we  have  been 
describing.  Furthermore,  this  system,  which  has  come  to  be  called  the 
"Lorenz  attractor"  is  closely  related  to  fluid  dynamic  equations  for 
convective  motions  in  a  heated  fluid.  Lorenz  chose  this  system  because 
it  is  simple  to  solve  and  related  to  more  general  fluid  dynamic  systems, 
such  as  those  described  by  the  Naviei — Stokes  equations. 

It  seems  to  be  generally  presumed  that  the  more  complicated  systems 
(and  nature  itself)  behave  in  a  fashion  which  is  at  least  qualitatively 
similar  to  Lorenz’  (1963)  simple  system,  which  was  later  extended  to  a 
somewhat  more  realistic  set  of  ordinary  differential  equations  derived 
from  the  shallow-water  equations  with  bottom  topography  (Lorenz,  1980). 
In  that  latter  paper,  Lorenz  concludes  with  some  remarks  about  the 
importance  of  attractors  and  strange  attractors  to  fluid  dynamics. 
Lorenz  states, 

"The  attractor  set  of  a  dynamical  system  is  for  practical  purposes 
the  set  of  points  in  phase  space  which  will  continue  to  be  encoun¬ 
tered  by  an  arbitrary  orbit  after  an  arbitrarily  long  time  has 
passed.  For  a  large  class  of  forced  dissipative  systems  the 
attractor  has  zero  volume,  i.e.,  an  arbitrarily  selected  point  in 
phase  space  is  almost  always  not  on  the  attractor.  When  the 
general  solution  is  aperiodic,  the  attractor  is  strange." 


Lorenz  (1980)  then  goes  on  to  discuss  the  degree  to  which  the  attractors 
of  simpler  equation  sets  approximate  the  strange  attractors  of  a  com¬ 
plete  set.  He  speculates  that  the  attractor  of  a  global  circulation 
model  with  100,000  variables  would  be  strange,  but  would  have  many  fewer 
dimensions  than  the  entire  phase  space,  perhaps  as  few  as  several  hun¬ 
dred  dimensions. 

The  connection  between  fractals  and  strange  attractors  takes  place 
in  phase  space.  Mandelbrot  (1983)  contends  "that  for  most  purposes  an 
attractor  is  strange  when  it  is  a  fractal."  According  to  Frisch  (1980) 
one  shortcoming  of  the  strange  attractor  approaches  is  that  they  seem 
capable  of  describing  only  temporal  chaos,  but  not  the  commonly  observed 
chaotic  spatial  structure  especially  at  small  scales.  He  suggests  that 
this  shortcoming  may  arise  from  the  phase  limited  spaces  used  in  most 
studies,  which  represent  the  system  on  a  coarse  grid,  or  with  limited 
spatial  Fourier  modes.  It  appears  likely  at  this  time  that  a  focus  on 
fractals  in  geometric  space,  rather  the  fractal  nature  of  the  strange 
attractors  of  fluid  systems  in  phase  space,  will  prove  more  fruitful  for 
practical  applications. 


5.  Observational  Evidence 


Procaccia  (1984)  gives  a  review  of  how  fractal  structures  and 
turbulence  might  affect  turbulent  diffusion,  fluctuations  of  passive 
scalars,  electromagnetic  wave  propagation  and  the  perimeters  of  clouds 
in  the  atmosphere.  He  examines  two  particle  (relative)  turbulent  diffu- 

sion.  He  considers  two  particles  (at  r1  and  r 2)  and  the  behavior  of  the 
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separation  distance  R  =  ^-r.,  Producecl  bY  their  relative  velocity  i.e. 
by  V  (t )  =  V.|  -  V2  : 

t 

R(t)  =  R  ( 0 )  +  /  V(x)  dx  (111-38) 

o 


for  isotropic  turbulence  <V  (t)>  =  0  so  that  the  ensemble  average  sepa¬ 
ration  (represented  by  angle  brackets,  <  >)  is  constant  and  equal  to  the 
initial  separation,  but  the  variance,  <R^>  changes: 


d_ 

dt 


<r2> 


t 

=2  J  <V(t)*V(x)>  dx 
o 


(II 1—39 ) 


The  above  equation  requires  estimation  of  a  time  correlation  function 
for  velocity  difference  over  a  length  scale  R.  Procaccia  (1984)  asserts 

that  although  the  correlation  <V(t)*V(x)>  is  nonstationary,  there  is 
some  function  of  scaled  time  variables  g(x)  such  that 

<V( t )  • V( x )>  =  <V(t)*V(t)>  g[(t-x )/t„] 


( I I 1-40 ) 


A 


K' 


where  tR  is  the  typical  decay  time  for  velocity  differences  over  a 
length  scale  R.  Substituting  in  Equation  (II 1—39 )  provides  asymptotic 
predictions 


It  <r2> 


<v(t)*v(t)>t 

<V(t)*V(t)>t 


U<tR 

t>>t 

R 


(Ill-ill ) 


at  extreme  times.  When  R  is  in  the  inertial  range,  the  diffusivity 
d<R2>/dt  can  be  determined  from  <V(t)*V(t)>  ,  which  is  relatively  easy 
to  determine.  For  the  "homogeneous  fractal  model"  of  turbulence 
described  earlier, 

<V(t)*V(t)>  -  <e>2/3  R2/3  (R/l  )1-D/3  (III-42) 

o 

where  X,  is  the  outer  scale  of  the  turbulence  and  tR  is  taken  to  be  the 
ratio  o?  the  separation  distance  R  to  a  typical  velocity  difference  over 
that  separation  distance.  Hentschel  and  Procaccia  (1983)  used  the  above 
assumption  to  derive  the  following  expressions  for  d<R2>/dt 


It  <r2> 


dt 


<R2> 


<£> 


1/3 


<e> 


1/3 


R,/3 

R*'3 


(R/fc) 

(R/l  ) 
o 


1/2-D/6 


2-2 D/3 


t<<t 

r 

t>>t 

r 


(III-43a) 

(III-U3b) 


Where  R  is  the  root-mean-square  of  the  separation.  For  space  filling 
turbulence,  where  D=3,  the  above  expressions  reduce  to  the  classical 
"4/3  law". 

Hentschel  and  Procaccia  (1983a)  examined  data  on  2-point  diffusion 
published  by  Richardson  (1926)  and  Gifford  (1957).  Their  results  from 
the  Gifford  data  give  an  estimate  for  D  between  2.5  and  2.75.  The 
Richardson  data  give  estimates  that  correspond  to  a  value  of  D  between 
2.64  and  2.78.  Hentschel  and  Procaccia  (1983)  concluded  that  these 
values  were  reasonable  and  supported  tne  "fractally  homogeneous  turbu¬ 
lence"  concept.  It  should  be  noted  that  they  did  not  exhaust  the  avail¬ 
able  data  and  that  further  tests  are  possible.  For  example,  Gifford 
(1977)  identified  21  sources  of  relative  diffusion  data. 

Procaccia  (1984)  also  discusses  Lovejoy's  studies  of  the  fractal 
properties  of  perimeters  of  clouds  and  rain  areas.  Lovejoy  (1982)  found 
that  the  perimeter  P  of  a  cloud  is  related  to  the  cloud  area,  A,  by 

P  -  (a)^/2  (III— 44) 


i 
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where  D  =  1.35  ±0.05.  Procaccia  argues  that  D  is  the  fractal  dimension 
of  the  perimeter,  and  that  the  fractal  dimension  of  the  cloud  surface 
Dc  =  1+D,  or  2.35  ±0.05,  (in  the  isotropic  case). 

Lovejoy's  (1982)  results  represented  a  range  of  six  orders  of 
magnitude  in  cloud  area  —  from  about  1  kilometer  to  1000  kilometers  in 
length.  Hentschel  and  Procaccia  (1984),  argue  that  the  perimeter  of  the 
cloud  is  on  a  surface  where  some  scalar  is  constant.  Lovejoy  (1982) 
defines  the  boundary  according  to  its  temperature.  Procaccia  (1984) 
considers  a  surface  defining  the  outer  boundary  of  the  cloud,  and  how 
that  surface  is  distorted  by  turbulent  diffusion.  He  then  relates  the 
fractal  dimension  of  the  resulting  surface  to  that  associated  with  rate 
of  change  of  separation  variance,  d<R2>/dt._  Expressed  in  terms  of  the 
fractal  dimension  of  the  cloud  perimeter,  6  the  relationship  is  (from 
Procaccia,  1984: 

5  =  ( 1 1 -D  )/6  .  (III-45) 

where  D  is  the  fractal  dimension  of  the  turbulence.  Using  the  values 
for  D,  that  Hentschel  and  Procaccia  (1983)  derived  from  Gifford's  (1957) 
data  gives  a  value  of  5  of  _about  1.4  which  is  certainly  in  reasonable 
agreement  with  the  value  of  6  =  1.35  ±0.05. 

From  the  standpoint  of  fluid  flow  modeling,  one  of  the  more  inter¬ 
esting  aspects  of  Hentschel  and  Procaccia’s  (1983,  1984;  Procaccia, 

1984)  analysis  is  the  link  it  provides  between  scalar  distributions  and 
the  properties  of  the  turbulent  field.  They  made  the  connection  via 
cloud  perimeters,  but  there  are  other  possibilities  that  deserve  explor¬ 
ation.  For  example,  the  distribution  of  smokes,  dyes  or  metal  flakes  in 
flow  visualization  experiments  might  be  analyzed  for  fractal  dimension 
and  used  to  estimate  turbulent  properties.  Ludwig  (1986)  has  been 
attempting  to  estimate  fractal  dimension  from  lidar  observations  of 
backscatter  in  vertical  cross  sections  through  smoke  plumes.  No  attempt 
has  been  made  to  determine  turbulent  properties,  but  that  seems  to  be 
the  logical  next  step. 


D.  Recent  Developments  in  Atmospheric  Applications 

To  this  point,  the  discussion  has  centered  on  isotropic  turbulence 
and  its  effects.  However,  horizontal  scales  of  atmospheric  motion 
extend  to  larger  dimensions  than  do  vertical  scales,  which  suggests  that 
any  derivations  that  depend  on  isotropic  turbulence  should  not  be  used 
throughout  the  atmosphere,  at  least  without  modification.  Recent  work 
by  Lovejoy  and  Schertzer  (1986)  and  Schertzer  and  Lovejoy  (1985)  has 
provided  the  basis  for  necessary  modification.  They  note  the  following 
important  facts  regarding  mesoscale  processes: 

_  g 

(1)  The  energy  spectrum  is  scaling  (i.e.  it  is  of  the  form  k  h 
where  k  is  wave  number  and  is  the  appropriate  value  for 

wave  numbers  in  the  horizontal  plane,  Bh  -  5/3). 
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(2)  The  energy  spectrum  for  wave  numbers  in  the  vertical  plane  is 
also  scaling,  but  anisotropy  makes  the  relevant  exponent,  Bv, 
quite  different;  -  11/5. 

(3)  There  is  extreme  variability,  because  active  regions  that 
account  for  most  of  the  energy  and  moisture  flux  are  very 
sparsely  distributed. 

Lovejoy  and  Schertzer  (1986)  present  most  of  the  concepts  related 
to  isotropic  turbulence  that  have  already  been  discussed  and  then  extend 
those  concepts  to  the  anisotropic  case.  They  make  the  extension  by 
continuing  to  assume  that  there  is  a  constant  energy  flux  over  the  range 
of  scales  of  interest,  and  that  there  are  rules  for  describing  how  the 
statistical  properties  of  eddies  are  transformed  from  one  scale  to 
another.  In  dealing  with  anisotropic  eddies,  their  shape  must  be 
characterized  along  with  their  dimensions.  Schertzer  and  Lovejoy  choose 
to  represent  the  anisotropic  cascade  in  a  manner  very  similar  to  that 
discussed  earlier  for  the  isotropic  case,  except  that  horizontal  and 
vertical  scaling  are  considered  separately,  i.e. 


Pr[Ac(AAx)  >  q]  =  Pr[A  Ac(Ax)  >  q] 


( III-1*  5a ) 


Pr [ Ac ( A  Ay )  >  q]  =  Pr[A  Ac ( Ay )  >  q] 


(III-1*  5b) 


Pr [ Ac ( A Az )  >  q]  »  Pr[A  Ac(Az)  >  q] 


(III-1*  5c) 


Equations  (III-1^)  can  be  written  in  matrix  form: 


Pr[ Ac (TAr )  >  q]  =  Pr[A  h  Ac ( Ar )  >  q]  , 


(III-46) 


The  matrix  T  is  given  by: 


A  0  0 


T  =  0  A  0 


(III-1*?) 


where 


0  0  A 


H  =  H  /H. 


(III—48) 


w 


(111-49) 


The  matrix  T  produces  a  magnification  overall,  with  stretching  in 
the  z  direction;  it  transforms  the  probability  distributions  and  intro¬ 
duces  an  elliptical  geometry  to  account  for  the  different  horizontal  and 
vertical  scalings.  Figure  D-6  shows  the  how  the  magnification  and 
stretching  relates  the  small,  vertically  oriented  eddies  to  the  a  large, 
horizontally  oriented  eddy.  The  transformation  changes  the  volume  of 
the  eddy  by  a  factor  A2  AHz  =  ADel,  where  the  "elliptical  dimension" 


Del  =  2  +  Hz-  In  an  isotropic  atmosphere,  Hz  =  1 ,  and  therefore  the 
elliptical  dimension  equals  3.  For  2-dimensionally  isotropic  turbulence 
Hz  =  0  and  Del  =  2.  By  analogy  with  the  purely  homogeneous  case, 
Schertzer  and  Love joy  (1985)  regard  Del  as  a  fractal  dimension  for  the 
anisotropic  space.  Schertzer  and  Lovejoy  (1985)  show  that  the  number  of 
eddies  of  horizontal  scale  l  is  proportional  to  H^el. 


The  reason  for  introducing  the  modification  described  above  is 
apparent  in  Figure  D-6;  it  provides  a  smooth  transition  from  the  very 
large,  horizontally  oriented,  "Hadley-like"  cells  down  to  the  vertically 
oriented  convective  cells.  There  is  some  scale  where  the  eddies  are — at 
least  in  the  statistical  sense--spherical .  Lovejoy  and  Schertzer  (1985) 
refer  to  this  as  the  "sphero-scale."  They  contend  that  it  is  quite 
variable  because  it  depends  on  both  the  average  turbulent  energy  flux 
e  and  the  average  flux  of  buoyant  force  variance. 


The  dependence  on  buoyant  force  flux  variance  arises,  according  to 
Schertzer  and  Lovejoy  (1985)  because  the  vertical  structure  is  governed 
by  that  factor,  analogous  to  the  flux  of  kinetic  energy  e.  The  flux  of 
the  buoyant  force  variance  D  is  given  by 


<|>(Az)  =  t  1  (Az)  Af2  ( Az) 


(II 1—50 ) 


where  t(Az)  is  a  characteristic  time  for  the  transfer  process,  and  f  is 
the  buoyant  force  per  unit  mass.  It  follows  from  dimensional  analysis 
that: 


Pr[Av(Az)  >  q]  =  Pr { [ $( Az ) j1 /5  Az3/5  >  q) 


( 1 11-51  ) 


According  to  Schertzer  and  Lovejoy  (1985),  the  above  scaling  holds  in 
the  vertical,  while  the  Kolmogorov  scaling  governs  the  horizontal,  i.e. 


Pr[Av(Ax)  >  q]  =  Pr  { [e  ( Ax  )  j1 /3  Ax1/3  >  q) 


(III-52 ) 


The  probability  distributions  for  the  two  Av's  can  be  equated  (after 
cubing)  to  give 


Pr(eAx  >  q)  =  PrU3^  Az^/^  >  q) 


(III-53) 
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Love joy  and  Schertzer  (1986) 


FIGURE  D-6  SCHEMATIC  DIAGRAM  ILLUSTRATING  LOVEJOY  AND  SCHERTZER' 
CONCEPT  OF  THE  CASCADE  FROM  LARGE  TO  SMALL  SCALES 


to  provide  an  algebraic  relationship  between  the  buoyant  force  variance 

and  eddy  kinetic  energy  fluxes.  For  the  case  with  hyperbolic  upper 

tails  of  the  probability  distributions,  and  exponents  a  and  a  for  the 

turbulent  energy  and  buoyant  force  variance  fluxes  respectively, 

Schertzer  and  Love joy  (1985)  show  that  a  -  5a, /3  . 

e  9 

One  of  the  more  important  consequences  of  Schertzer  and  Lovejoy’s 
(1985)  analysis  is  that  it  provides  a  bridge  between  the  large  and  small 
scales  of  motion.  It  does  not  require  a  distinction  between  large  scale 
motions  which  are  essentially  2-dimensional,  and  small  motions  which  are 
more  3-dimensional.  They  contend  that  there  is  no  evidence  (although  it 
has  been  sought)  for  any  transition  from  2-dimensional  to  3-dimensional 
regimes.  It  should  be  noted  that  their  hypothesis  is  not  necessarily 
relevant  only  to  the  atmosphere,  but  may  also  be  relevant  to  other  fluid 
systems  where  there  are  vastly  different  horizontal  and  vertical  scales, 
and  where  vertical  motions  are  largely  the  product  of  buoyancy  effects. 


IV  APPLICATION  TO  CLOSURE  SCHEMES  AND  FLOW  MODELING 


It  is  beyond  the  scope  of  this  review  to  develop  a  turbulent  clo¬ 
sure  scheme  for  flow  modeling  based  on  fractal  concepts.  The  literature 
does  contain  some  tantalizing  results  that  appear  to  have  the  potential 
to  solve  some  of  closure  problems.  However,  the  whole  history  of  turbu¬ 
lence  modeling  is  characterized  by  promising  approaches  that  have  just 
one  essential  item  beyond  grasp.  In  essence,  the  problem  is  to  deter¬ 
mine  local  covariances  between  important  hydrodynamic  variables  (e.g. 
velocity  components,  temperature,  pressure,  diffusing  scalar  concentra¬ 
tions  and  so  forth)  for  volumes  that  are  representati ve  the  model’s  grid 
dimensions.  Determination  of  these  covariances  must  not  require  infor¬ 
mation  that  is  not  generated  by  the  model  for  the  larger  scales  (i.e.  it 
must  use  data  available  at  the  grid  points). 

The  first  thing  that  is  apparent  from  a  review  of  the  fractal 
literature  is  that  it  should  be  possible  to  generate  scalar  and  vector 
fields  that  have  realistic  statistical  properties.  However,  this  does 
not  by  itself  solve  the  closure  problem.  Two  more  major  requirements 
remain: 

(1)  it  must  be  possible  to  infer  the  fractal  characteristics  of 
those  smaller  scale  distributions  from  the  information  at 
model  grid  points. 

(2)  The  determination  of  covariances  between  different  quantities 
(especially  anisotropic  covariance)  requires  both  that  the 
statistics  of  the  variable  distributions  in  space  be  realis¬ 
tic,  and  that  their  mutual  distributions  be  correctly  related. 


Is  there  anything  in  the  literature  that  suggests  that  the  above 
problems  can  be  resolved?  Until  the  last  two  or  three  years,  the  answer 
would  have  been  no.  However,  the  contention  by  Schertzer  and  Lovejoy 
(1985)  that  there  is  a  continuous  cascade  from  the  synoptic  atmospheric 
scales  down  to  the  viscous  scale,  and  that  the  statistical  properties  of 
that  cascade  can  be  described  using  the  concept  of  elliptical  fractal 
dimension  suggests  that  it  should  be  possible  effect  a  transition  from 
the  grid  scale  to  smaller  scales.  Furthermore,  the  relationships  among 
elliptical  fractal  dimension,  turbulent  energy  flux  and  the  flux  of 
buoyant  force  variance  involve  major  hydrodynamic  variables,  which  gives 
hope  that  any  subgrid  scale  fields  that  might  be  generated  would  have 
the  appropriate  interrelationships  among  velocity,  temperature,  pres¬ 
sure,  and  so  forth. 

The  work  of  Hentschel  and  Procaccia  (1983)  attempts  to  relate  the 
fractal  properties  of  diffusing  scalars  to  the  fractal  dimension  of  the 
turbulence.  That  work  might  also  provide  the  link  necessary  to  obtain 
the  covariances.  Certainly,  it  has  the  potential  to  help  in  the  design 
of  useful  experimental  studies,  because  it  should  considerably  easier  to 
measure  small  scale  distributions  of  a  scalar  in  space  than  to  measure 
the  distributions  of  velocity  fluctuations. 

The  intermittency  of  turbulence  means  that  on  all  scales  there  are 
"active"  and  "inactive"  regions.  It  seems  reasonable  to  presume  that 
subgrid  scale  effects  in  regions  that  are  active  on  scales  resolved  by 
the  model  will  be  different  from  those  in  areas  which  are  inactive.  We 
might  also  expect  that  the  processes  occuring  in  active  regions  will 
smooth  out  fluctuations,  so  that  those  regions  become  less  active.  The 
smaller  velocity  fluctuations  characteristic  of  the  inactive  regions  may 
allow  the  development  of  gradients  that  will  in  turn  lead  to  greater 
fluctuations  and  cause  the  region  to  eventually  become  "active." 
Regardless  of  whether  or  not  this  conjecture  is  true,  it  is  clear  that 
any  model  should  be  able  to  deal  with  time  variations  of  intermittency 
in  the  turbulent  fields. 

Two  major  components  will  be  required  for  any  closure  scheme 
derived  from  fractal  concepts.  They  are: 

(1  )  A  method  for  deriving  fractal  characteristics  and  identifying 
"active"  areas  from  information  at  the  model  gridpoints. 

(2)  A  method  for  determining  the  necessary  covariance  terms  from 
the  fractal  characteristics,  either  analytically  or  through 
parameterization  of  numerical  "experiments"  using  artificially 
generated  fields  based  on  specified  fractal  characteristics. 

There  are  major  unsolved  problems  associated  with  both  the  above 
components.  Determining  fractal  characteristics  requires  calculations 
of  spatial  correlations  over  a  wide  range  of  scales.  Furthermore, 
anisotropy  requires  that  the  correlations  in  the  vertical  and  horizontal 
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planes  be  determined  separately.  The  limited  number  of  gridpoints 
available  in  most  modeling  may  make  this  unreliable  or  impossible.  This 
may  not  be  as  severe  a  problem  as  it  first  appears,  because  the  litera¬ 
ture  provides,  through  theory  or  observation,  estimated  values  for  many 
of  the  scaling  exponents  and  probability  distribution  parameters.  It  is 
encouraging  to  note  that  the  parameters  of  interest  in  turbulence  model¬ 
ing  such  as  turbulent  kinetic  energy,  buoyant  flux  variance,  velocity 
fluctuations,  and  so  forth  have  generally  been  the  focus  of  the  studies 
reviewed  here.  Thus,  the  most  severe  problems  may  be  limited  to  the 
determination  of  the  degree  of  turbulent  activity  in  an  area  which  could 
be  a  much  simpler  task. 

One  final  problem  remains.  Most  of  the  derivations  and  analyses  do 
not  consider  boundaries  like  the  earth's  surface  or  the  basin  and  air- 
water  interfaces  of  water  bodies.  This  may  limit  the  usefulness  of  much 
published  information,  at  least  in  the  vicinity  of  the  boundaries. 


V  CONCLUDING  REMARKS 


It  is  clear  that  a  great  deal  of  effort  will  be  required  to  exploit 
whatever  potential  there  is  for  application  of  fractal  concepts  to  the 
modeling  of  turbulence.  Other  fractal  applications  are  much  easier  and 
have  in  fact  been  accomplished.  Visual  simulation  of  cloud  and  smoke 
plume  appearance  has  been  done  (e.g.  Medler  et  al,  1986;  Love  joy  and 
Schertzer,  1986).  It  is  a  fairly  easy  step  from  there  to  the  determina¬ 
tion  of  concentration  probability  distributions,  effects  on  optical 
transmission  and  so  forth.  Of  course,  it  requires  that  the  fractal 
characteristics  of  the  temperature  field,  the  diffusing  smoke  cloud  or 
whatever  scalar  is  of  interest  be  known,  but  this  should  be  experiment¬ 
ally  determinable.  In  the  case  of  atmospheric  cloud  formations  and 
precipitation  areas  the  fractal  properties  have  been  determined.  The 
information  is  already  at  hand  for  solving  many  practical  problems  such 
as  those  related  to  scintillation  and  other  optical  transmission 
effects,  and  those  related  to  variations  in  the  visibility  through 
fog. 
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The  fact  that  Schertzer  and  Love joy  (1986)  have  made  connections 
between  buoyant  flux  variance  and  turbulent  kinetic  energy  lends  cre¬ 
dence  to  the  notion  that  one  can  use  fractal  concepts  to  model  correla¬ 
tion  functions  as  well  as  variance  of  a  single  variable.  It  is  of 
course  the  correlation  functions  which  are  the  key  to  the  modeling  of 
turbulent  effects  in  fluid  flow. 

Even  the  problem  of  boundaries  may  be  tractable.  According  to 
Mandelbrot  (1983)  natural  surfaces  assume  fractal  shapes.  It,  there¬ 
fore,  seems  plausible  that  air  motions  above  such  a  surface  would  assume 
fractal  properties.  There  remains  of  course  the  problem  of  transition 
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between  those  motions  where  the  surface  plays  an  important  role  and 
those  aloft  were  the  variations  are  introduced  by  turbulent  kinetic 
energy  cascades  and  buoyancy  forces.  Even  if  fractal  concepts  could 
only  provide  turbulence  models  in  the  free  fluid,  away  from  boundaries, 
that  would  be  a  major  contribution. 
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